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Abstract: The method described here for recovering the shape of a surface from a 
shaded image can deal with complex, wrinkled surfaces. Integrability can be enforced 
easily because both surface height and gradient are represented (A gradient field is 
integrable if it is the gradient of some surface height function). The robustness of the 
method stems in part from linearization of the reflectance map about the current 
estimate of the surface orientation at each picture cell (The reflectance map gives the 
dependence of scene radiance on the gradient; or, equivalently, surface orientation). 
The new scheme can find an exact solution of a given shape-from- shading problem 
even though a regularizing term is included. The reason is that the penalty term is 
needed only to stabilize the iterative scheme when it is far from the correct solution; 
it can be turned off as the solution is approached. This is a reflection of the fact 
that shape-from-shading problems are not ill-posed when boundary conditions are 
available or when the image contains singular points. 

This report includes a review of previous work on shape from shading and pho- 
toclinometry. Novel features of the new scheme are introduced one at a time to make 
it easier to see what each contributes. Included is a discussion of implementation de- 
tails that are important if exact algebraic solutions of synthetic shape-from-shading 
problems are to be obtained. The hope is that better performance on synthetic data 
will lead to better performance on real data. 
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1. Background 

The first method developed for solving a shape-from-shading problem was 
restricted to surfaces with special reflecting properties [Rindfleisch 66]. For the 
surfaces that Rindfleisch considered, profiles of the solution can be obtained 
by integrating along predetermined straight lines in the image. The general 
problem was formulated and solved later [Horn 70, 75], using the method 
of characteristic strip expansion applied to the nonlinear first-order partial 
differential image irradiance equation. When the light sources and the viewer 
are far away from the scene being viewed, use of the reflectance map makes 
the analysis of shape-from-shading algorithms much easier [Horn 77] [Horn 
& Sjoberg 79]. Several iterative schemes, mostly based on minimization of 
some functional containing an integral of the brightness error, arose later 
[Woodham 77] [Strat 79] [Ikeuchi & Horn 81] [Kirk 84, 87] [Brooks & Horn 85] 
[Horn & Brooks 86] [Frankot & Chellappa 88]. 

The new method presented here was developed in part as a response 
to recent attention to integrability 1 [Horn k Brooks 86] [Frankot & Chel- 
lappa 88] and exploits the idea of a coupled system of equations for depth and 
slope [Harris 86, 87] [Horn 88]. It borrows from well-known variational ap- 
proaches to the problem [Ikeuchi & Horn 81] [Brooks & Horn 85] and an exist- 
ing least-squares method for estimating surface shape given a needle map (see 
[Ikeuchi 84], chapter 11 in [Horn 86], and [Horn & Brooks 86]). For one choice 
of parameters, the new method becomes similar to one of the first iterative 
methods ever developed for shape from shading on a regular grid [Strat 79], 
while it degenerates into another well-known method [Ikeuchi & Horn 81] for a 
different choice of parameters. If the brightness error term is dropped, then it 
becomes a surface interpolation method [Harris 86, 87]. The computational ef- 
fort grows rapidly with image size, so the new method can benefit from proper 
multigrid implementation [Brandt 77, 80, 84] [Brandt & Dinar 79], as can 
existing iterative shape-from-shading schemes [Terzopolous 83, 84] [Kirk 84, 
87]. Alternatively, one can apply so-called direct methods for solving Poisson's 
equations [Simchony, Chellappa & Shao 89]. 

Linear expansion of the reflectance map about the current estimate of 
the surface gradient leads to more rapid convergence. More importantly, this 
modification often allows the scheme to converge when simpler schemes di- 
verge, or get stuck in local minima of the functional. Most existing iterative 
shape-from-shading methods handle only relatively simple surfaces and so 
could benefit from a retrofit of this idea. 

The new scheme was tested on a number of synthetic images of increas- 
ing complexity, including some generated from digital terrain models of steep, 

*A gradient field is integrable if it is the gradient of some surface height function. 
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Figure 1. Reconstruction of surface from shaded image. See text. 
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wrinkled surfaces, such as a glacial cirque with numerous gullies. Shown in 
Figure 1(a) is a shaded view of a digital terrain model, with lighting from 
the Northwest. This is the input provided to the algorithm. The underlying 
231 x 178 digital terrain model was constructed from a detailed contour map, 
shown in Figure 2, of Huntington ravine on the eastern slopes of Mount Wash- 
ington in the White Mountains of New Hampshire 2 . Shown in Figure 1(b) 
is a shaded view of the same digital terrain model with lighting from the 
Northeast. This is not available to the algorithm, but is shown here to make 
apparent features of the surface that may not stand out as well in the other 
shaded view. Figure 1(c) shows a shaded view of the surface reconstructed 
by the algorithm, with lighting from the Northwest— it matches Figure 1(a) 
exactly. More importantly, the shaded of view of the reconstructed surface 
with lighting from the Northeast, shown in Figure 1(d), matches Figure 1(b) 

exactly 3 . 

With proper boundary conditions, the new scheme recovers surface ori- 
entation exactly when presented with noise-free synthetic scenes. Previous 
iterative schemes do not find the exact solution, and in fact wander away 
from the correct solution when it is used as the initial guess. To obtain exact 
algebraic solutions, several details of the implementation have to be care- 
fully thought through, as discussed in section 6. Simple surfaces are easier to 
process — with good results even when several of the implementation choices 
are not made in an optimal way. Similarly, these details perhaps may be of 
lesser importance for real images, where other error sources could dominate. 

In the next few sections we review image formation and other elementary 
ideas underlying the usual formulation of the shape-from-shading problem. 
Photoclinometry is also briefly reviewed for the benefit of researchers in ma- 
chine vision who may not be familiar with this field. We then discuss both 
the original and the variational approach to the shape-from-shading problem. 
Readers familiar with the basic concepts may wish to skip over this material 
and go directly to section 5, where the new scheme is derived. For additional 
details see chapters 10 and 11 in Robot Vision [Horn 86] and the collection of 
papers, Shape from Shading [Horn &; Brooks 89]. 

2. Review of Problem Formulation 

2.1 Image Projection and Image Irradiance 

For many problems in machine vision it is convenient to use a camera-centered 
coordinate system with the origin at the center of projection and the Z-axis 

2 The gullies are steep enough to be of interest to ice-climbers. 

3 For additional examples of reconstructions from shaded images, see section 7. 
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Figure 2. Contour map from which the digital terrain model used 
to synthesize Figures 1(a) and (b) was interpolated. The surface was 
modelled as a thin plate constrained to pass through the contours at the 
specified elevations. The interpolating surface was found by solving the 
biharmonic equation, as described at the end of section 5.4. 



aligned with the optical axis (the perpendicular from the center of projection 
to the image plane) 4 . We can align the X- and F-axes with the image plane 
x- and y-axes. Let the principal distance (that is, the perpendicular distance 
from the center of projection to the image plane) be /, and let the image plane 
be reflected through the center of projection so that we avoid sign reversal of 



4 In photoclinometry it is customary to use an object- centered coordinate system. 
This is because surface shape can be computed along profiles only when strong 
additional constraint is provided, and such constraints are best expressed in an 
object-centered coordinate system. Working in an object-centered coordinate 
system, however, makes the formulation of the shape-from-shading problem con- 
siderably more complex (see, for example, [Rindfleisch 66]). 
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the coordinates. Then the perspective projection equations are 

x = /— and y = f~z- C 1 ) 

The problem is simplified if we assume that the depth range is small compared 
with the distance of the scene from the viewer (which is often the case when 
we have a narrow field of view, that is, when we use a telephoto lens). Then 
we have 

xnl-X and y«^F, (2) 

for some constant Zo, so that the projection is approximately orthographic. 
In this case it is convenient to rescale the image coordinates so that we can 
write x = X and y = Y. For work on shape from shading it is also convenient 
to use z, height above some reference plane perpendicular to the optical axis, 
rather than the distance measured along the optical axis from the center of 
projection. 

If we ignore vignetting and other imaging system defects, then image irra- 
diance E at the point (x, y) is related to scene radiance L at the corresponding 
point in the scene by [Horn 86] 

£ = L^(j) 2 cos 4 a, (3) 

where d is the diameter of the lens aperture, / is the principal distance, and 
the off-axis angle a is given by 

tan a = — yjx 2 + y 2 . (4) 

Accordingly, image irradiance 5 is a multiple of the scene radiance, with the 
factor of proportionality depending inversely on the square of the /-number 6 . 
If we have a narrow field of view, the dependence on the off-axis angle a 
can be neglected. Alternatively, we can normalize the image by dividing 
the observed image irradiance by cos 4 a (or whatever the actual vignetting 
function happens to be). 

We conclude from the above that what we measure in the image is directly 
proportional to scene radiance, which in turn depends on the strength and 
distribution of illumination sources, the surface micro-structure and surface 
orientation. 

In order to be able to solve the shape from shading problem from a single 
image we must assume that the surface is uniform in its reflecting properties. 
If we also assume that the light sources are far away, then the irradiance of 



5 Grey-levels are quantized estimates of image irradiance. 

6 The /-number is the ratio of the principal distance to the diameter of the aperture, 
that is, f/d. 
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different parts of the scene will be approximately the same and the incident 
direction may be taken as constant. Finally, if we assume that the viewer 
is far away, then the direction to the viewer will be roughly the same for all 
points in the scene. Given the above, we find that scene radiance does not 
depend on the position in space of a surface patch, but only on its orientation. 

2.2 Specifying Surface Orientation 

Methods for recovering shape from shading depend on assumptions about 
the continuity of surface height and its partial derivatives. First of all, since 
shading depends only on surface orientation, we must assume that the surface 
is continuous and that its first partial derivatives exist. Most formulations 
implicitly also require that the first partial derivatives be continuous, and 
some even require that second partial derivatives exist. The existence and 
continuity of derivatives lends a certain "smoothness" to the surface and allows 
us to construct local tangent planes. We can then talk about the local surface 
orientation in terms of the orientation of these tangent planes. 

There are several commonly used ways of specifying the orientation of a 
planar surface patch, including: 

• Unit surface normal n [Horn & Brooks 86]; 

• Point on the Gaussian sphere [Horn 84]; 

• Surface gradient (p, q) [Horn 77]; 

• Stereographic coordinates (/,#) [Ikeuchi & Horn 81]; 

• Dip and strike (as defined in geology) 7 ; 

• Luminance longitude and latitude 8 ; 

• Incident and emittance angles (i and e) 9 ; 

For our purposes here, the components of the surface gradient 

p = Tx and q = % (5) 



7 Dip is the angle between a given surface and the horizontal plane, while strike is 
the direction of the intersection of the surface and the horizontal plane. The line 
of intersection is perpendicular to the direction of steepest descent. 

8 Luminance longitude and latitude are the longitude and latitude of a point on 
a sphere with the given orientation, measured in a spherical coordinate system 
with the poles at right angles to both the direction towards the source and the 
direction toward the viewer. 

9 Incident and emittance angles are meaningful quantities only when there is a single 
source; and even then there is a two-way ambiguity in surface orientation unless 
additional information is provided. The same applies to luminance longitude and 
latitude. 
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will be most directly useful for specifying surface orientation. 

We can convert between different representations easily. For example, 
suppose that we are to determine the unit surface normal given the gradient 
components. We know that if we move a small distance 8x in x, then the 
change in height is 8z — pSx (since p is the slope of the surface in the x 
direction). Thus (1,0, p) T is a tangent to the surface. If we move a small 
distance Sy in y, then the change in height is 8z = q 8y (since q is the slope 
of the surface in the y direction). Thus (0,1, q) T is also a tangent to the 
surface. The normal is perpendicular to all tangents, thus parallel to the 
cross-product of these particular tangents, that is parallel to (— p, —q, 1) . 
Hence a unit normal can be written in the form 

ft= ; u ',, 2 (-*>--^> T - <«> 

y/l+p 2 + q 2 
Note that this assumes that the z-component of the surface normal is positive. 
This is not really a problem since we can only see surface elements whose 
normal vectors point within 7r/2 of the direction toward the viewer — other 
surface elements are turned away from the viewer. 

We can use the same notation to specify the direction to a collimated 
light source or a small portion of an extended source. We simply give the 
orientation of a surface element that lies perpendicular to the incident light 
rays. So we can write 10 

for some p a and q a . 

2.3 Reflectance Map 

We can show the dependence of scene radiance on surface orientation in the 
form of the reflectance map R(p, q). The reflectance map can be depicted 
graphically in gradient space 11 as a series of nested contours of constant 
brightness [Horn 77, 86]. 

The reflectance map may be determined experimentally by mounting a 
sample of the surface on a goniometer stage and measuring its brightness un- 
der the given illuminating conditions for various orientations. Alternatively, 



10 There is a small problem, however, with this method for specifying the direction 
toward the light source: A source may be "behind" the scene, with the direction 
to the source more than n/2 away from the direction toward the viewer. In this 
case the z-component of the vector pointing toward the light source is negative. 

11 The coordinates of gradient space are p and q, the slopes of the surface in the x 
and y direction respectively. 
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one may use the image of a calibration object (such as a sphere) for which 
surface orientation is easily calculated at every point. Finally, a reflectance 
map may be derived from a phenomenological model, such as that of a Lam- 
bertian surface. In this case one can integrate the product of the bidirectional 
reflectance distribution function (BRDF) and the given distribution of source 
brightness as a function of incident angle [Horn &; Sjoberg 79]. 

An ideal Lambertian surface illuminated by a single point source provides 
a convenient example of a reflectance map 12 . Here the scene radiance is given 
by R(p, q) = (E /k) cos «', where i is the incident angle (the angle between the 
surface normal and the direction toward the source), while E is the irradiance 
from the source on a surface oriented perpendicular to the incident rays. (The 
above formula only applies when i < 7r/2; the scene radiance is, of course, zero 
for i > 7r/2.) Now cos i = n • s, so 

T>( \ E * 1+PsP + qsQ /ON 

7T y/1 + p 2 + q 2 y/l + p 2 s + q 2 s 
as long as the numerator is positive, otherwise R(p, q) = 0. 



2.4 Image Irradiance Equation 

We are now ready to write down the image irradiance equation 

E(x,y) = f3R(p(x,y),q(x,y)), (9) 

where E(x,y) is the irradiance at the point (x,y) in the image, while R(p,q) 
is the radiance at the corresponding point in the scene, at which p = p(x,y) 
and q = q(x, y). The proportionality factor fl depends on the /-number of the 
imaging system (and may include a scaling factor that depends on the units 
in which the instrument measures brightness). It is customary to rescale 
image irradiance so that this proportionality factor may be dropped. If the 
reflectance map has a unique global extremum, for example, then the image 
can be normalized in this fashion, provided that a point can be located that 
has the corresponding surface orientation 13 . 

Scene radiance also depends on the irradiance of the scene and a re- 
flectance factor (loosely called albedo here). These factors of proportionality 



12 Note that shape-from-shading methods are most definitely not restricted to Lam- 
bertian surfaces. Such special surfaces merely provide a convenient pedagogical 
device for illustrating basic concepts. 

13 If there is a unique maximum in reflected brightness, it is convenient to rescale 
the measurements so that this extremum corresponds to E = 1. The same ap- 
plies when there is a unique minimum, as is the case for the scanning electron 
microscope (SEM). 
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can be combined into one that can be taken care of by normalization of im- 
age brightness. Then only the geometric dependence of image brightness on 
surface orientation remains in R(p, q), and we can write the image irradiance 
equation in the simple form 

E{x,y) = R{p{x,y),q{x,y)) (10) 

or 

E(x,y) = R(z x (x,y),z 9 (x,y)), (11) 

where p — z x and q = z y are the first partial derivatives of z with respect 
to x and y. This is a first-order partial differential equation that is typically 
nonlinear, because the reflectance map in most cases depends nonlinearly on 
the gradient. 

2.5 Reflectance Map Linear in Gradient 

Viewed from a sufficiently great distance, the material in the maria of the 
moon has the interesting property that its brightness depends only on lu- 
minance longitude, being independent of luminance latitude [Hapke 63, 65]. 
When luminance longitude and latitude are related to the incident and emit- 
tance angles, it is found that longitude is a function of (cos i/ cos e). From the 
above we see that cosi = fi • s, while cose = n • v, where v = (0,0, l) r is a 
unit vector in the direction toward the viewer. Consequently, 

cosi n-s 1 , s. / no x 

= ^r^r = —, - A l+PsP + q s q). (12) 

cos e n • v y/\ + p\ + q\ 

Thus (cosi/ cose) depends linearly on the gradient components p and q, and 
we can write 

R(p,q) = f(cp + sq), (13) 

for some function / and coefficients c and s. Both Lommel-Seeliger's and 
Hapke's functions fit this mold [Minnaert 61] [Hapke 63, 65]. (For a few other 
papers on the reflecting properties of surfaces, see [Hapke 81, 84] [Hapke & 
Wells 81] and the bibliography in [Horn k Brooks 89].) We can, without loss 
of generality 14 , arrange for c 2 + s 2 = 1. 

If the function / is continuous and monotonic 15 , we can find an inverse 

cp + sq = f- 1 (E(x,y)). (14) 



14 We see that c : s = p s : q a , so that the direction specified in the image by (c,s) is 
the direction "toward the source," that is, the projection into the image plane of 
the vector s toward the light source. 

15 If the function / is not monotonic, there will be more than one solution for 
certain brightness values. In this case one may need to introduce assumptions 
about continuity of the derivatives in order to decide which solution to choose. 
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The slope in the image direction (c, s) is 

m = i^i = - 7 =L=/-'(B(x,y)). (15) 

Vc 2 + s 2 vc 2 + s 2 

We can integrate 16 out this slope along the line 

x(0 = x + c£ and y(£) = y + s £, (16) 

to obtain 

z(0 = z + -=L= /* /" * (j5(*0?), yfo))) d»y. (17) 

Vc 2 + s 2 ./o v ' 

An extension of the above approach allows one to take into account perspective 

projection as well as finite distance to the light source [Rindfleisch 66]. Two 

changes need to be made; one is that the reflectance map now is no longer 

independent of image position (since the directions to the viewer and the 

source vary significantly); and the other that the integral is for the logarithm 

of the radial distance from the center of projection, as opposed to distance 

measured parallel to the optical axis. 

The above was the first shape-from-shading or photoclinometric problem 
ever solved in other than a heuristic fashion. The original formulation was 
considerably more complex than described above, as the result of the use of full 
perspective projection, the lack of the notion of anything like the reflectance 
map, and the use of an object-centered coordinate system [Rindfleisch 66]. 

Note that we obtain profiles of the surface by integrating along predeter- 
mined straight lines in the image. Each profile has its own unknown constant 
of integration, so there is a great deal of ambiguity in the recovery of surface 
shape. In fact, if z(x,y) is a solution, so is 

z(x,y) = z(x,y)+g(sx-cy) (18) 

for an arbitrary function g\ This is true because 

z x = z x +sg'(sx-cy) and z y = z y - cg'(sx - cy), (19) 

so 

cp + sq = cp + sq, (20) 

where p = z x and q - z y . It follows that R(p, q) = R(p, q). This ambiguity can 
be removed if an initial curve is given from which the profiles can be started. 
Such an initial curve is typically not available in practice. Ambiguity is not 
restricted to the special case of a reflectance map that is linear in the gradient: 
Without additional constraint shape-from-shading problems typically do not 
have a unique solution. 



16 The integration is, of course, carried out numerically, since the integrand is derived 
from image measurements and not represented as an analytic function. 
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2.6 Low Gradient Terrain and Oblique Illumination 

If we are looking at a surface where the gradient (p, q) is small, we can ap- 
proximate the reflectance map using series expansion 

R(p, q) « J2(0, 0) + p J2p(0, 0) + q Rg(0, 0). (21) 

This approach does not work when the reflectance map is rotationally sym- 
metric, since the first-order terms then drop out 17 . If the illumination is 
oblique, however, we can apply the method in the previous section to get a 
first estimate of the surface. Letting c = R p (0,0), s = i2 g (0, 0) and 

f- 1 (E(x,y))=E(x,y)-R(0,0), (22) 

we find that 

z(0 = z + 1 / (^(zfo), yfo)) - R(0, 0)) dq. (23) 

^(0,0) + Rl(0,0)Jo V 

(For a related frequency domain approach see [Pentland 88].) 

One might imagine that the above would provide a good way to get initial 
conditions for an iterative shape from shading method. Unfortunately, this is 
not very helpful, because of the remaining ambiguity in the direction at right 
angles to that of profile integration. Iterative methods already rapidly get 
adequate variations in height along "down-sun profiles," but then struggle for 
a long time to try to get these profiles tied together in the direction at right 
angles. 

The above also suggests that errors in gradients of a computed solution 
are likely to be small in the direction towards or away "from the source" 
and large in the direction at right angles. It should also be clear that it is 
relatively easy to find solutions for slowly undulating surfaces (where p and q 
remain small) with oblique illumination (as in [Kirk 87]). It is harder to deal 
with cases where the surface gradient varies widely, and with cases where the 
source is near the viewer. 

3. Brief Review of Photoclinometry 

Photoclinometry is the recovery of surface slopes from images [Wilhelms 64] 
[Rindfleisch 66] [Lambiotte & Taylor 67] [Watson 68] [Lucchitta & Gambell 70] 
[Tyler, Simpson & Moore 71] [Rowan, McCauley & Holm 71] [Bonner & 
Small 73] [Wildey 75] [Squyres 81] [Howard, Blasius & Cutt 82]. Many papers 



17 The reflectance map is rotationally symmetric, for example, when the source is 
where the viewer is, or when an extended source is symmetrically distributed 
about the direction toward the viewer. 
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and abstracts relating to this subject appear in places that may seem inacces- 
sible to someone working in machine vision [Davis, Soderblom, & Eliason 82] 
[Passey & Shoemaker 82] [Davis & McEwen 84] [Davis & Soderblom 83, 84] 
[Malin & Danielson 84] [Wilson et al. 84] [McEwen 85] [Wilson et al. 85] (For 
additional references see [Horn k Brooks 89]). Superficially, photoclinometry 
may appear to be just another name for shape from shading. Two different 
groups of researchers independently tackled the problem of recovering surface 
shape from spatial brightness variations in single images. Astrogeologists and 
workers in machine vision became aware of each other's interests only a few 
years ago. The underlying goals of the two groups are related, but there are 
some differences in approach that may be worthy of a brief discussion. 

3.1 Photoclinometry versus Shape from Shading 

• First, photoclinometry has focused mostly on profile methods (photo- 
clinometrists now refer to existing shape-from-shading methods as area- 
based photoclinometry, as opposed to -profile-based). This came about 
in large part because several of the surfaces of interest to the astroge- 
ologist have reflecting properties that allow numerical integration along 
predetermined lines in the image, as discussed above in section 2.5 [Rind- 
fieisch 66]. Later, a similar profile integration approach was applied to 
other kinds of surfaces, by using strong assumptions about local surface 
geometry instead. The assumption that the surface is locally cylindrical 
leads to such a profile integration scheme [Wildey 86], for example. More 
commonly, however, it has been assumed that the cross-track slope is 
zero, in a suitable object-centered coordinate system [Squyres 81]. This 
may be reasonable when one is considering a cross-section of a linearly ex- 
tended feature, like a ridge, a graben, or a central section of a rotationally 
symmetric feature like a crater. 

• The introduction of constraints that are easiest to express in an object- 
centered coordinate system leads away from use of a camera-centered 
coordinate system and to complex coordinate transformations that tend 
to obscure the underlying problem. A classic paper on photoclinometry 
[Rindfleisch 66] is difficult to read for this reason, and as a result had little 
impact on the field. On the other hand, it must be acknowledged that 
this paper dealt properly with perspective projection, which is important 
when the field of view is large. In all but the earliest work on shape 
from shading [Horn 70, 75], the assumption is made that the projection 
is approximately orthographic. This simplifies the equations and allows 
introduction of the reflectance map. 
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• The inherent ambiguity of the problem does not stand out as obviously 
when one works with profiles, as it does when one tries to fully reconstruct 
surfaces. This is perhaps why workers on shape from shading have been 
more concerned with ambiguity, and why they have emphasized the im- 
portance of singular points and occluding boundaries [Bruss 82] [Deift & 
Sylvester 81] [Brooks 83] [Blake, Zisserman k Knowles 85] [Saxberg 88]. 

• The recovery of shape is more complex than the computation of a set of 
profiles. Consequently much of the work in shape from shading has been 
restricted to simple shapes. At the same time, there has been extensive 
testing of shape from shading algorithms on synthetic data. This is some- 
thing that is important for work on shape from shading, but makes little 
sense for the study of simple profile methods, except to test for errors in 
the procedures used for inverting the photometric function. 

• Shape-from-shading methods easily deal with arbitrary collections of col- 
limated light sources and extended sources, since these can be accommo- 
dated in the reflectance map by integrating the BRDF and the source 
distribution. In astrogeology there is only one source of light (if we ig- 
nore mutual illumination or interflection between surfaces), so methods 
for dealing with multiple sources or extended sources were not developed. 

• Calibration objects are used both in photoclinometry and shape from 
shading. In photoclinometry the data derived is used to fit parameters to 
phenomenological models such as those of Minnaert, Lommel and Seel- 
iger, Hapke, and Lambert. In work on shape from shading the numerical 
data is at times used directly without further curve fitting. The pa- 
rameterized models have the advantage that they permit extrapolation 
of observations to situations not encountered on the calibration object. 
This is not an issue if the calibration object contains surface elements 
with all possible orientations, as it will if it is smooth and convex. 

• Normalization of brightness measurements is treated slightly differently 
too. If the imaging device is linear, one is looking for a single overall 
scale factor. In photoclinometry this factor is often estimated by looking 
for a region that is more or less flat and has known orientation in the 
object-centered coordinate system. In shape from shading the brightness 
of singular points is often used to normalize brightness measurements 
instead. The choice depends in part on what is known about the scene, 
what the shapes of the objects are (that is, are singular points or oc- 
cluding boundaries imaged) and how the surface reflects light (that is, is 
there a unique global extremum in brightness). 

• Finally, simple profiling methods usually only require continuity of the 
surface and existence of the first derivative (unless there is an ambiguity 



14 Height and Gradient from Shading 

in the inversion of the photometric function whose resolution requires that 
neighboring derivatives are similar). Most shape from shading methods 
require continuous first derivatives and the existence of second derivatives 
(In some cases use is made of the equality of the second cross-derivatives 
taken in different order, that is, z xy = z yx ). This means that these meth- 
ods do not work well on scenes composed of objects that are only piecewise 
smooth, unless appropriately modified 18 (but see [Malik & Maydan 89]). 

3.2 Profiling Methods 

We have seen in section 2.5 how special photometric properties sometimes 
allow one to calculate a profile by integration along predetermined straight 
lines in the image. The other approach commonly used in photoclinometry 
to permit simple integration is to make strong assumptions about the surface 
shape, most commonly that, in a suitable object-centered coordinate system, 
the slope of the surface is zero in a direction at right angles to the direction in 
which the profile is being computed. Local surface orientation has two degrees 
of freedom. The measured brightness provides one constraint. A second 
constraint is needed to obtain a solution for surface orientation. A known 
tangent of the surface can provide the needed information. Two common 
cases are treated in astrogeology: 

(a) features that appear to be linearly extended (such as some ridges and 
grabens), in a direction presumed to be "horizontal" (that is, in the local 
average tangent plane); 

(b) features that appear to be rotationally symmetric (like craters), with 
symmetry axis presumed to be "vertical" (that is, perpendicular to the 
average local tangent plane). 

In each case, the profile is taken "across" the feature, that is, in a direction 
perpendicular to the intersection of the surface with the average local tangent 
plane. Equivalently, it is assumed that the cross-track slope is zero in the 
object-centered coordinate system. 

One problem with this approach is that we obtain a profile in a plane 
containing the viewer and the light source, not a "vertical" profile, one that is 
perpendicular to the average local tangent plane. One way to deal with this 
is to iteratively adjust for the image displacement resulting from fluctuations 
in height on the surface, using first a scan that really is just a straight line in 



18 Methods for recovering the shapes of polyhedral objects using shading on the faces 
and the directions of the projections of the edges into the image are discussed in 
[Sugihara 86] and [Horn 86]. 
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the image, then using the estimated profile to introduce appropriate lateral 
displacements into the scan line, and so on [Davis & Soderblom 84]. 

It turns out that the standard photoclinometric profile approach can 
be easily generalized to arbitrary tangent directions, ones that need not be 
perpendicular to the profile, and also to nonzero slopes. All that we need 
to assume is that the surface can locally be approximated by a (general) 
cylinder, that is, a surface generated by sweeping a line, the generator, along 
a curve in space. Suppose the direction of the generator is given by the vector 
t = (a, b,c) T . Note that at each point on the surface, a line parallel to the 
generator is tangent to the surface. Then, since the normal is perpendicular 
to any tangent, we have t • n = at every point on the surface, or just 

ap + bq = c. (24) 

This, together with the equation E = R(p, q), constitutes a pair of equations 
in the two unknowns p and q. There may, however, be more than one solution 
(or perhaps none) since one of the equations is nonlinear. Other means must 
be found to remove possible ambiguity arising from this circumstance. Under 
appropriate oblique lighting conditions, there will usually only be one solution 
for most observed brightness values. 

From the above we conclude that we can recover surface orientation lo- 
cally if we assume that the surface is cylindrical, with known direction of the 
generator. We can integrate out the resulting gradient in any direction we 
please, not necessarily across the feature. Also, the generator need not lie 
in the local tangent plane; we can deal with other situations, as long as we 
know the direction of the generator in the camera-centered coordinate system. 
Further generalizations are possible, since any means of providing one more 
constraint on p and q will do. 

In machine vision too, some workers have used strong local assumptions 
about the surface to allow direct recovery of surface orientation. For example, 
if the surface is assumed to be locally spherical, the first two partial derivatives 
of brightness allow one to recover the surface orientation [Pentland 84] [Lee 
& Rosenfeld 85]. Alternatively, one may assume that the surface is locally 
cylindrical [Wildey 84, 86] to resolve the ambiguity present locally in the 
general case. 

4. Review of Shape from Shading Schemes 

4.1 Characteristic Strips 

The original solution of the general shape from shading problem uses the 
method of characteristic strip expansion [Horn 70, 75]. The basic idea is quite 
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easy to explain using the reflectance map [Horn 77, 86]. Suppose that we are 
at a point (x, y, z) T on the surface and we wish to extend the solution a small 
distance in some direction by taking a step 8x in x and 8y in y. We need to 
compute the change in height 8z. This we can do if we know the components 
of the gradient, because 

8z=p8x + q8y. (25) 

So, as we explore the surface, we need to keep track of p and q in addition to 
x, y and z. This means that we also need to be able to compute the changes 
in p and q when we take the step. This can be done using 

8p = r 8x + s8y and 8q = s 8x + 1 8y, (26) 

where r = z xx , s = z xy = z yx and t — z yy are the second partial derivatives of 
the height. It seems that we need to now keep track of the second derivatives 
also, and in order to do that we need the third partial derivatives, and so on. 

To avoid this infinite recurrence, we take another tack. Note that we 
have not yet used the image irradiance equation E(x,y) = R(p,q). To find 
the brightness gradient we differentiate this equation with respect to x and y 
and so obtain 

E x =rR p + sR q and E y =sR p +tR q . (27) 

At this point we exploit the fact that we are free to choose the direction of 
the step (8x,8y). Suppose that we pick 

8x = R p 8£ and 8y = R q 8£, (28) 

then, from equations (26) & (27) we have 

8p = E x 8i and 8q = E y 8(. (29) 

This is the whole "trick." We can summarize the above in the set of ordinary 
differential equations 

x = R p , y = R q , z—pR p + qR q 

p = E x , q = E y , (30) 

where the dot denotes differentiation with respect to £, a parameter that varies 
along a particular solution curve (the equations can be rescaled to make this 
parameter be arc length). Note that we actually have more than a mere 
characteristic curve, since we also know the orientation of the surface at all 
points in this curve. This is why a particular solution is called a characteristic 
strip. The projection of a characteristic curve into the image plane is called 
a base characteristic. 

The base characteristics are predetermined straight lines in the image 
only when the ratio x : y = R p : R q is fixed, that is when the reflectance 
map is linear in p and q. In general, one cannot integrate along arbitrary 
curves in the image. Also, an initial curve is needed from which to sprout the 
characteristics strips. 
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It turns out that direct numerical implementation of the above equations 
does not yield particularly good results, since the paths of the characteristics 
are affected by noise in the image brightness measurements and errors tend 
to accumulate along their length. In particularly bad cases, the base charac- 
teristics may even cross, which does not make any sense in terms of surface 
shape. It is possible, however, to grow characteristic strips in parallel and use 
a so-called sharpening process to keep neighboring characteristics consistent 
[Horn 70, 75]. This greatly improves the accuracy of the solution, since the 
computation of surface orientation is tied more closely to image brightness 
itself rather than to the brightness gradient. This also makes it possible to 
interpolate new characteristic strips when existing ones spread too far apart, 
and to remove some when they approach each other too closely. 

4.2 Rotationally Symmetric Reflectance Maps 

One can get some idea of how the characteristics explore a surface by con- 
sidering the special case of a rotationally reflectance map, as might apply 
when the light source is at the viewer (or when dealing with scanning electron 
microscope (SEM) images). Suppose that 

R(p,q) = f(p 2 + q 2 ), (3i) 

then 

R P = 2p /V + q 2 ) and R q = 2q f'(p 2 + q 2 ), (32) 

and so the directions in which the base characteristics grow are given by 

x = k p and y = kq, (33) 

for some k. That is, in this case the characteristics are curves of steepest 
ascent or descent on the surface. The extrema of surface height are sources and 
sinks of characteristics. These are the points where the surface has maxima 
in brightness. 

This example illustrates the importance of so-called singular points. At 
most image points, as we have seen, the gradient is not fully constrained by 
image brightness. Now suppose that R(p, q) has a unique global maximum 19 , 
that is 

R(p, q) < R(p , q ) for all (p, q) ± (p , q ) . (34) 

A singular point (zo, yo) in the image is a point where 

E(x ,yo) = R(po,qo)- (35) 

At such a point we may conclude that (p,q) = (po,9o)- Singular points in 
general are sources and sinks of characteristic curves. Singular points provide 



19 The same argument applies when the unique extremum is a minimum, as it is in 
the case of scanning electron microscope (SEM) images. 
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strong constraint on possible solutions [Horn 70, 75] [Bruss 82] [Brooks 83] 
[Saxberg 88]. 

A limb is the image of an occluding boundary, where the surface is tan- 
gent to the direction toward the viewer. It has been suggested that occluding 
boundaries provide strong constraint on possible solutions [Ikeuchi & Horn 81] 
[Bruss 82]. As a consequence there has been interest in representations for 
surface orientation that behave well near the occluding boundary, unlike the 
gradient, which becomes infinite [Ikeuchi & Horn 81] [Horn & Brooks 86]. 
Recently there has been some question as to how much constraint occlud- 
ing boundaries really provide, given that singular points appear to already 
strongly constrain the solution [Brooks 83] [Saxberg 88]. 

4.3 Existence and Uniqueness 

Amazingly, questions of existence and uniqueness of solutions of the shape- 
from-shading problem have still not been resolved satisfactorily. One problem 
is that it is hard to say anything in general that applies to all possible re- 
flectance maps. More can be said when specific reflectance maps are chosen, 
such as ones that are linear in the gradient or those that are rotationally 
symmetric [Bruss 82]. 

It has recently been shown that there exist impossible shaded images, 
that is, images that do not correspond to any surface illuminated in the spec- 
ified way [Szeliski & Horn 89]. It may turn out that almost all images with 
multiple singular points are impossible in this sense [Saxberg 88]. This is an 
important issue, because it may help explain how our visual system sometimes 
determines that the surface being viewed cannot possibly be uniform in its 
reflecting properties. One can easily come up with smoothly shaded images, 
for example, that do not yield an impression of shape, instead appearing as 
flat surfaces with spatially varying reflectance or surface "albedo." (See also 
Figure 10 in section 7.2). 

4.4 Variational Formulations 

As discussed above in section 2.4, in the case of a surface with constant albedo, 
when both the observer and the light sources are far away, surface radiance 
depends only on surface orientation and not on position in space and the 
image projection can be considered to be orthographic 20 . In this case the 



20 The shape-from-shading problem can be formulated and solved when the viewer 
and the light sources are not at a great distance [Rindfleisch 66] [Horn 70, 75], 
but then scene radiance depends on position as well as surface orientation, and 
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image irradiance equation becomes just 

E(x, y) = R(p(x, y), q(x, y)) , (36) 

where E(x,y) is the image irradiance at the point (x,y), while R(p,q), the 
reflectance map, is the (normalized) scene radiance of a surface patch with 
orientation specified by the partial derivatives 

p=|i and g = |, (37) 

of surface height z(x,y) above some reference plane perpendicular to the op- 
tical axis. 

The task is to find z(x,y) given the image E(x,y) and the reflectance 
map R{p, q). Additional constraints, such as boundary conditions and sin- 
gular points, are needed to ensure that there is a unique solution [Bruss 82] 
[Deift & Sylvester 81] [Blake, Zisserman & Knowles 85] [Saxberg 88]. If we 
ignore integrability 21 , some versions of the problem of shape from shading 
may be considered to be ill-posed 22 , that is, there is not a unique solution 
{p(x,y),q(x,y)} that minimizes the brightness error 

(E(x,y)-R(p,q)) 2 dxdy. (38) 

In fact the error can be made equal to zero for an infinite number of choices for 
{p(x,y), q(x, y)}. We can choose one of these solutions by finding the one that 
minimizes some functional such as a measure of "departure from smoothness" 

JJ(pl+P 2 y + Ql+^)dxdy, (39) 

while satisfying the constraint E(x,y) = R(p,q). Introducing a Lagrange 
multiplier \(x,y) to enforce the constraint, we find that we have to minimize 



//< 
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((Pi + P\ + ql + <& + A(x, y) (E - R)) dx dy. (40) 

The Euler equations are 

Ap+X(x,y)R p =0 and Aq + \(x,y)R q = 0. (41) 

After elimination of the Lagrange multiplier A(x,y), we are left with the pair 
of equations 

R q Ap = R p Aq and E(x,y) = R(p, q). (42) 

There may not be any simple convergent iterative scheme for this constrained 
variational problem [Horn & Brooks 86] (compare [Wildey 75]). 



the notion of a reflectance map is not directly applicable. 
21 A gradient-field (or needle-diagram) {p(x,y),q(x,y)} is integrable if there exists 

some surface height function z(x,y) such that p(x,y) = z x (x,y) and q(x,y) = 

z y (x,y), where the subscripts denote partial derivatives. 
22 If we impose integrability, and provide suitable boundary conditions, then the 

shape-from-shading problem is definitely not ill-posed. 
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7/o 



We can approach this problem another way using the "departure from 
smoothness" measure in a penalty term [Ikeuchi & Horn 81], looking instead 
for a minimum of 

'(E(x,y) - R(p,q)) 2 + X(p 2 x + p 2 y + q 2 x + ??)) dxdy. (43) 

It should be pointed out that a solution of this "regularized" problem is not a 
solution of the original problem, although it may be close to some solution of 
the original problem [Brooks 85]. In any case, this variational problem leads 
to the following coupled pair of second-order partial differential equations: 

A Ap = - (E(x, y) - R(p, q)) R p (p, q) 

XAq = -(E(x,y)-R(p,q))R q (p,q) 
Using a discrete approximation of the Laplacian operator 23 

{A/}*, « £(J„ - /„), (45) 

where / is a local average 24 of /, and e is the spacing between picture cells, 
we arrive at the set of equations 

k X' pki = K X' p kl + (E(x , y ) - R(p, q)) R p (p, q) 

nX'qu = KX'q kl + (E(x,y) - R(p, q)) R q (p, q) 

where A' = A/e 2 . This set of equations suggests the iterative scheme 

Pil +1) =Al ) + -\j(E(x,y) - R(p^,q^))R p (p^,q^) 

kX> V ) (47) 

?& +1) = i$ + ^(E(x,y) - R( p ( n \q^))R q (p^\q^) 

where the superscript denotes the iteration number 25 . 

From the above it may appear that R(p, <?), R p (p, q), and R g (p, q) should 
be evaluated using the "old" values of p and q. It turns out that the numerical 
stability of the scheme is somewhat enhanced if they are evaluated instead at 
the local average values p and q [Ikeuchi & Horn 81]. 

One might hope that the correct solution of the original shape-from- 
shading problem provides a stable point for the iterative scheme. This is not 
too likely, however, since we are solving a modified problem that includes a 



There are several methods for approximating the Laplacian operator, including 
five-point and nine-point approximations. It is well known that, while the nine- 
point approximation involves more computation, its lowest-order error term has 
a higher order than that of the five-point approximation. 
24 Here k = 4 when the local average f k i is computed using the four edge-adjacent 
neighbors, while k = 10/3, when 1/5 of the average of the corner-adjacent neigh- 
bors is added to 4/5 of the average of the edge-adjacent neighbors (see later). 
These equations are solved iteratively because the system of equations is so large 
and because of the fact that the reflectance map R(p,q) is typically nonlinear. 
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penalty term. Consequently, an interesting question one might ask about an 
algorithm such as this, is whether it will "walk away" from the correct solution 
of the original image irradiance equation E(x,y) = R(p, q) when this solution 
is provided as an initial condition [Brooks 85]. The algorithm described here 
does just that, since it can trade off a small amount of brightness error against 
an increase in surface smoothness. At the solution, we have E(x, y) = R(p, q), 
so that the right hand sides of the two coupled partial differential equations 
(equations (44)) are zero. This implies that if the solution of the modified 
problem is to be equal the solution of the original problem then the Laplacians 
of p and q must be equal to zero. This is the case for very few surfaces, just 
those for which 

Az(x,y) = k, (48) 

for some constant k. While this includes all harmonic functions, it excludes 
most real surfaces, for which adjustments away from the correct shape are 
needed to assure equality of the left and right hand sides of equations (44) 
describing the solution of the modified problem. In general, this approach pro- 
duces solutions that are too smooth, with the amount of distortion depending 
on the choice of the parameter A. For related reasons, this algorithm does 
well only on simple smooth shapes, and does not perform well on complex, 
wrinkled surfaces. 

4.5 Recovering Height from Gradient 

In any case, we are also still faced with the problem of dealing with the lack 
of integrability, that is the lack of a surface z(x,y) such that p(x,y) — z x (x, y) 
and q(x,y) = z y (x,y) 26 . We can try to find the surface that has partial 
derivatives z x and z y that come closest to matching the computed p(x,y) and 
q(x,y) by minimizing 

jj((z x -p) 2 + (z y -q) 2 )dxdy. (49) 

This leads to the Poisson equation 

Az=p x + q y . (50) 

Using the discrete approximation of the Laplacian given above (equation (45)) 
yields 

-=■*« = -rZfcf - (Px + q y ), (51) 

a set of equations that suggests the following iterative scheme: 



„(»+i) _ ~(«) 



z kl = z 



ff " ^ ({p,}&> + {*,}$) , (52) 



26 The resulting gradient field is likely not to be integrable because we have not 
enforced the condition p y = q x , which corresponds to z xy = z yx . 
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where the terms in braces are numerical estimates of the indicated derivatives 
at the picture cell (fc, /). 

The so-called natural boundary conditions 27 here are just 

cz x +sz y = cp + sq, (53) 

where (c, s) is a normal to the boundary. 

Another way of dealing with the integrability issue is to directly minimize 

'{E{x,y) - R{p,q)) 2 + \(p y - q x f) dxdy. (54) 

This leads to the coupled partial differential equations [Horn & Brooks 86] 

*(Pyy ~ Qxy) = (E(x,y) - R(p,q))R p , 

Mtf*z -Pyx) = (E(x,y) - R(p,q))R q . 
This set of equations can also be discretized by introducing appropriate finite 
difference approximations for the second partial derivatives p yy , q xx and the 
cross derivatives of p and q. An iterative scheme is suggested once one isolates 
the center terms of the discrete approximations of p yy and q xx . This is very 
similar to the method developed by Strat, although he arrived at his scheme 
directly in the discrete domain [Strat 79]. His iterative scheme avoids the 
excessive smoothing of the one described early, but appears to be less stable, 
in the sense that it diverges under a wider set of circumstances. 

5. New Coupled Height and Gradient Scheme 

The new shape-from- shading scheme will be presented through a series of 
increasingly more robust variational methods. We start with the simplest, 
which grows naturally out of what was discussed in the previous section. 

5.1 Fusing Height and Gradient Recovery 

One way of fusing the recovery of gradient from shading with the recovery of 
height from gradient, is to represent both gradient (p, q) and height z in one 
variational scheme and to minimize the functional 



jJ({E{x,y) - R{p,q)) 2 + fi((z x -pf + (z y - qf)) dxdy. 



(56) 



Note that, as far as p(x,y) and q(x, y), are concerned this is an ordinary calcu- 
lus problem (since no partial derivatives of p and q appear in the integrand). 
Differentiating the integrand with respect to p(x,y) and q(x,y) and setting 



Natural boundary conditions arise in variational problems where no boundary 
conditions are explicitly imposed [Courant Sz Hilbert 62]. 
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the result equal to zero leads to 

P = z x + -(E - R)R P , 

\ (57) 

q = z y + -{E - R)R q . 

Now z(x,y) does not occur directly in (E(x,y) — R(p,q)) so we actually just 
need to minimize 

JJ((z x - pf + (z, - q) 2 ) dx dy, (58) 

and we know from the previous section that the Euler equation for this vari- 
ational problem is just 

Az =p x + q y . (59) 

We now have one equation for each of p, q and z. 

These three equations are clearly satisfied when p = z x , q = z y and 
E — R. That is, if a solution of the original shape-from-shading problem 
exists, then it satisfies this system of equations exactly (which is more than 
can be said for some other systems of equations obtained using a variational 
approach, as pointed out in section 4.4). It is instructive to substitute the 
expressions obtained for p and q in p x + q y '■ 

Pz + Qy = z« + z yy + —{(E - R)(R PP p x + R pq (p y + 9x) + R qq q y ) (60) 

- (R 2 pPx + R p R q (py + q x ) + R 2 q q y ) + (E X R P + E y R q )Bigr). 

Since Az = (p x + (fy), we note that the three equations above are satisfied 
when 

(R 2 pPx + R p R q (p y +q x ) + R 2 q q y ) - (E x R p + E y R q ) (61) 

= (E- R)(R pp p x + R pg (Py + q x ) + R qq q y ). 

This is exactly the equation obtained at the end of section 4.2 in [Horn & 
Brooks 86], where an attempt was made to directly impose integrability us- 
ing the constraint p y = q x . It was stated there that no convergent iterative 
scheme had been found for solving this complicated nonlinear partial differ- 
ential equation directly. The method described here provides an indirect way 
of solving this equation. 

Note that the natural boundary conditions for z are once again 

cz x + sz y = cp + sq, (62) 

where (c, s) is a normal to the boundary. 

The coupled system of equations above for p, q and z immediately sug- 
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gests an iterative scheme 

pfr" = Mtf + Ue - R)R,, 
V- 

qll^^Mff + hE-R)^, (63) 

4? +,) = 4" ) + 7({p.)t? +,) + te}i? +1) ). 

where we have used the discrete approximation of the Laplacian for z intro- 
duced in equation (45). This new iterative scheme works well when the initial 
values given for p, q and z are close to the solution. It will converge to the ex- 
act solution if it exists; that is, if there exist a discrete set of values {zjfc/} such 
that {pu} and {qu} are the discrete estimate of the first partial derivatives 
of z with respect to x and y respectively and Eki — R(pkl,qu)- 

In this case the functional we wish to minimize can actually be reduced 
to zero. It should be apparent that for this to happen, the estimator used for 
the Laplacian must match the sum of the convolution of the discrete estimator 
of the x derivative with itself and the convolution of the discrete estimator of 
the y derivative with itself. (This and related matters are taken up again in 
section 6.2.) 

The algorithm can easily be tested using synthetic height data zki- One 
merely estimates the partial derivatives using suitable discrete difference for- 
mulae and then uses the resulting values pki and qui to compute the synthetic 
image Eki = R(pkl,qkl)- This construction guarantees that there will be an 
exact solution. If a real image is used, there is no guarantee that there is an 
exact solution, and the algorithm can at best find a good discrete approx- 
imation of the solution of the underlying continuous problem. In this case 
the functional will in fact not be reduced exactly to zero. In some cases the 
residue may be large. This may be the result of aliasing introduced when 
sampling the image, as discussed in section 6.5, or because in fact the image 
given could not have arisen from shading on a homogeneous surface with the 
reflectance properties and lighting as encoded in the reflectance map — that 
is, it is an impossible shaded image. 

The iterative algorithm described in this section, while simple, is not very 
stable unless one is close to the exact solution, particularly when the surface 
is complex and the reflectance map is not close to linear in the gradient. It 
can be improved greatly by linearizing the reflectance map. It can also be 
stabilized by adding a penalty term for departure from smoothness. This 
allows one to come close to the correct solution, at which point the penalty 
term is removed in order to prevent it from distorting the solution. We first 
treat the linearization of the reflectance map. 



5. New Coupled Height and Gradient Scheme 25 

5.2 Linearization of Reflectance Map 

We can develop a better scheme than the one developed in the previous sec- 
tion, while preserving the apparent linearity of the equations, by approximat- 
ing the reflectance map R(p, q) locally by a linear function of p and q. There 
are several options for choice of reference gradient for the series expansion, so 
let us keep it general for now at (po,9o) 28 - We have 

R(p, q) « R(p , q ) + (p- Po) R P (po,qo) + (q - qo) R g (po,qo) + • • • (64) 
Again, gathering all of the terms in pu and qu on the left hand sides of the 
equations, we now obtain 

(fi + R 2 p )pki + RpRq qu = fiz x + (E-R- p R p - q R q )Rp, 
R q RpPki + (v + R 2 q )qki = nz y + (E -R- p Rp - q R q )R q , 
while the equation for z remains unchanged. (Note that here R, R p and R q 
denote quantities evaluated at the reference gradient (po,qo)). 

It is convenient to rewrite these equations in terms of quantities relative 
to the reference gradient (po,qo). Let 

Spu = Pki - Po and (fyfci = qu - qo, 

(66) 
Sz x = z x -po and Sz y = z y - q . K 

This leads to 

(fi + R 2 p ) Spu + RpR q Sq kl = fi6z x + (E- R)R P , 

R P Rq Sq k i + (n + R 2 q )Sq kl = fi6z y + (E- R)R q . 

(The equations clearly simplify somewhat if we choose (z x , z y ) for the reference 
gradient (po, qo).) We can view the above as a pair of linear equations for Spu 
and Sqki. The determinant of the 2x2 coefficient matrix, 

D = fi(fi + R 2 p + R 2 ), (68) 

is always positive, so there is no problem with singularities. The solution is 
given by 

DS P ki = (fi + R 2 )A- R p R q B, 

DSq kl = (fi + R 2 p )B-R q RpA, 

where 

A = fi8z x + (E-R)R P , 

B = fi6z y + (E-R)R q . 

This leads to a convenient iterative scheme where the new values are given by 

pi n , +1) =P ( o n) +^P k r l ) -d q ^ +1) = q ^+6qi!\ (71) 



The reference gradient will, of course, be different at every picture cell, but to 
avoid having subscripts on the subscripts, we will simple denote the reference 
gradient at a particular picture cell by (po,qo). 
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in terms of the old reference gradient and the increments computed above. 
The new version of the iterative scheme does not require a great deal more 
computation than the simpler scheme introduced in section 4.5, since the 
partial derivatives R p and R q are needed in any case. 

5.3 Incorporating Departure from Smoothness Term 

We now introduce a penalty term for departure from smoothness, effectively 
combining the iterative method of [Ikeuchi & Horn 81] for recovering p and 
q from E(x,y) and R(p,q), with the scheme for recovering z given p and q 
discussed in section 4.5. (For the moment we do not linearize the reflectance 
map; this will be addressed in section 5.6). We look directly for a minimum 
of 

"\E(x,y)-R(p,q)) 2 



//« 



+ Hpl+p 2 y + ql + q 2 9 ) (72) 

+ fi((z x -p) 2 + (z y -q) 2 )) dxdy. 

The Euler equations of this calculus of variations problem lead to the following 
coupled system of second-order partial differential equations: 

A Ap = -(E - R)R P - ft(z x - p), 

XAq = -(E - R)R q - fi(z s - q), (73) 

Az = p x + q y . 
A discrete approximation of these equations can be obtained by using the 
discrete approximation of the Laplacian operator introduced in equation (45): 

~Z (Pki ~ Pu) = ~{E - R)R P - n(z x - p k i), 



kX 

K 



(Iki ~ Ml) = ~(E - R)R q - K*y ~ Ml), ( 74 ) 



2 (zu - Zki) =p x +q y . 

where E, R, R p , and R q are the corresponding values at the picture cell (k, I), 

while z x , z y , p x and q y are discrete estimates of the partial derivative of z, p 

and q there. We can collect the terms in pu, qu and Zki on one side to obtain 

(kX' + n)p kl = (nX'p kl + fiz x ) + (E- R)R P , 

(kX' + ft) q kl = (kX' q kl + fiz y ) + (E- R)R q , ( 75 ) 

where A' = A/e 2 . These equations immediately suggest an iterative scheme, 
where the right hand sides are computed using the current values of the z k i, 



5. New Coupled Height and Gradient Scheme 27 

Pki, and qu, with the results then used to supply new values for the unknowns 
appearing on the left hand sides 

From the above it may appear that R(p, q), R p (p, <jf), and R q (p, q) should 
be evaluated using the "old" values of p and q. One might, on the other hand, 
argue that the local average values p and q, or perhaps even the gradient 
estimates z x and z y , are more appropriate. Experimentation suggests that 
the scheme is most stable when the local averages p and q are used. 

The above scheme contains a penalty term for departure from smooth- 
ness, so it may appear that it cannot to converge to the exact solution. Indeed, 
it appears as if the iterative scheme will "walk away" from the correct solution 
when it is presented with it as initial conditions, as discussed in section 4.4. It 
turns out, however, that the penalty term is needed only to prevent instability 
when far from the solution. When we come closer to the solution, A' can be 
reduced to zero, and so the penalty term drops out. It is tempting to leave 
the penalty term out right from the start, since this simplifies the equations 
a great deal, as shown in section 5.1. The contribution from the penalty term 
does, however, help damp out instabilities when far from the solution and so 
should be included. This is particularly important with real data, where one 
cannot expect to find an exact solution. 

(Note also that the coupled second order partial differential equations 
above (equation (76)) are eminently suited for solution by coupled resistive 
grids [Horn 88].) 

5.4 Relationship to Existing Techniques 

• Recently a new method has been developed that combines the itera- 
tive scheme discussed in section 4.4 for recovering surface orientation 
from shading with a projection onto the subspace of integrable gradients 
[Frankot & Chellappa 88]. The approach there is to alternately take one 
step of the iterative scheme [Ikeuchi k Horn 81] and to find the "nearest" 
integrable gradient. This gradient is then provided as initial conditions 
for the next step of the iterative scheme, ensuring that the gradient field 
never departs far from integrability. The integrable gradient closest to a 
given gradient field is found using orthonormal series expansion and by 
exploiting the fact that differentiation in the spatial domain corresponds 
to multiplication by frequency in the transform domain. 

• Similar results can be obtained by using instead the method described 
in section 4.5 for recovering the height z(x,y) that best matches a given 
gradient. The resulting surface can then be (numerically) differentiated 
to obtain initial values for p(a;,y) and q(x,y) for the next step of the 
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iterative scheme. This works, but not as well as the new scheme described 
in the previous section. 

• Next, note that we obtain the scheme of [Ikeuchi & Horn 81] (who ignore 
the integrability problem) discussed in section 4.4, if we drop the depar- 
ture from integrability term in the integrand — that is, when p = 0. If we 
instead remove the departure from smoothness term in the integrand — 
that is, when A = — we obtain something reminiscent of the iterative 
scheme of [Strat 79], although Strat dealt with the integrability issue in 
a different way. 

• Finally, if we drop the brightness error term in the integrand, we obtain 
the scheme of [Harris 86, 87] for interpolating from depth and slope. He 
minimizes 

// (A (Pi +P 2 y + ql + q 2 y ) + ((*, - P? + (z, -q) 2 )) dx dy. (76) 

and arrives at the Euler equations 

A Ap = -(z x - p), 

\Aq = -(z,-q), (77) 

Az = p x + q y . 

Now consider that 

A(Az) = A( Px + q y ). (78) 

Since application of the Laplacian operator and differentiation commute 
we have 

A(Az) = (A P ) x + (Aq) y , (79) 



or 



A A(Az) = -(z xx - Px ) - (z yy - q y ), (80) 



and so 



A A(Az) = -Az + ( Px + q y ) = 0. (81) 

So his method solves the biharmonic equation for z, by solving a coupled 
set of second-order partial differential equations. It does it in an elegant, 
stable way that permits introduction of constraints on both height z and 
gradient (p, q). This is a good method for interpolating from sparse depth 
and surface orientation data. 

The biharmonic equation has been employed to interpolate digital terrain 
models (DTMs) from contour maps. Such DTMs were used, for example, 
in [Horn &: Bachman 78] [Horn 81] [Sjoberg & Horn 83]. The obvious im- 
plementations of finite difference approximations of the biharmonic operator, 
however, tend to be unstable because some of the weights are negative, and 
because the corresponding coefficient matrix lacks diagonal dominance. Also, 
the treatment of boundary conditions is complicated by the fact that the 
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support of the biharmonic operator is so large. The scheme described above 
circumvents both of these difficulties — it was used to interpolate the digital 
terrain model used for the example illustrated by Figure 1. 

5.5 Boundary Conditions &; Nonlinearity of Reflectance Map 

So far we have assumed that suitable boundary conditions are available, that 
is, the gradient is known on the boundary of the image region to which the 
computation is to be applied. If this is not the case, the solution is likely not 
to be unique. We may nevertheless try to find a solution by imposing so-called 
natural boundary conditions [Courant & Hilbert 62]. The natural boundary 
conditions for the variational problem described here can be shown to be 

cp x + sp y = and c q x + s q y = (82) 

and 

cz x + sz y = cp + sq (83) 

where (c, s) is a normal to the boundary. That is, the normal derivative of 
the gradient is zero and the normal derivative of the height has to match the 
slope in the normal direction computed from the gradient. 

In the above we have approximated the original partial differential equa- 
tions by a set of discrete equations, three for every picture cell (one each for 
p, q and z). If these equations were linear, we could directly apply all the 
existing theory relating to convergence of various iterative schemes and how 
one solves such equations efficiently, given that the corresponding coefficient 
matrices are sparse 29 . Unfortunately, the equations are in general not linear, 
because of the nonlinear dependence of the reflectance map R(p, q) on the gra- 
dient. In fact, in deriving the above simple iterative scheme, we have treated 
R(p, q), and its derivatives, as constant (independent of p and q) during any 
particular iterative adjustment of p and q. 

5.6 Local Linear Approximation of Reflectance Map 

In section 5.2 we linearized the reflectance map in order to improve the simple 
scheme developed in section 5.1. We now do the same for the more complex 
scheme described in section 5.4. We have 

R(p, q) « R(p , q ) + (p- Po)R P (po,qo) + (q - ?o) R g (po,qo) + ••• (84) 
Again, gathering all of the terms in pki and qu on the left hand sides of the 
equations, we now obtain 

(X" + R 2 p )p k , + R p R q q kl 

29 See [Lee 88] for a proof of convergence of an iterative shape- from-shading scheme. 
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= (k\'pu + fiz x ) + (E - R- poR p - q Rq)R p , 
R q R p Pki + (\" + R 2 q )q k i (85) 

= (K\'q kl + fiZy) + (E - R — p Rp — q Rq)R q , 

while the equation for z remains unchanged. (Note that here R, R p and R q 
again denote quantities evaluated at the reference gradient (po,<7o))- In the 
above we have abbreviated A" = kA' + //. 

It is convenient to rewrite these equations in terms of quantities defined 
relative to the reference gradient: 

6pkl = Pki ~ Po and 8q k \ = qu - 9o 
S Pkl = Pkl ~ Po and 8q kl = q kl - q (86) 

Sz x = z x — po and 6z y = z y — qo 
This yields 

(A" + R\) 6 Pk , + R p R g Sq kl = k\' 6p kl + n6z x + (E- R)R P , 
RpR q 6q k i + (X" + R 2 q )6q k i^K\'6q kl + n6z y + (E-R)R q . 

(The equations clearly simplify somewhat if we choose either p and q or z x 
and z y for the reference gradient po and 50-) We can view the above as a pair 
of linear equations for 8p k i and Sq k \. The determinant of the 2x2 coefficient 
matrix 

D = \"(X" + R 2 p + R 2 q ) (88) 

is always positive, so there is no problem with singularities. The solution is 
given by 

DS Pkl = (\" + R 2 q )A-RpR q B, 

D6q kl = (\" + R 2 p )B-R q R p A, 

where 

A = k\' Sp kl + n6z x +(E- R)R P , 

B = kX' Sq kl + nSz y +(E- R)R q . 

This leads to a convenient iterative scheme where the new values are given by 

Pil +1) =pi n) +^kl ) and ^7 +1) = «{ ,l) + ^? ) , (91) 

in terms of the old reference gradient and the increments computed above. It 
has been determined empirically that this scheme converges under a far wider 
set of circumstances than the one presented in the previous section. 

Experimentation with different reference gradients, including the old val- 
ues of p and q, the local average p and q, as well as z x and z y showed that 
the accuracy of the solution and the convergence is affected by this choice. It 
became apparent that if we do not want the scheme to "walk away" from the 
correct solution, then we should use the old value of p and q for the reference 
po and q . 
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6. Some Implementation Details 

6.1 Derivative Estimators and Staggered Grids 



In one dimension, it is well-known from numerical analysis that the best finite 
difference estimators of even derivatives have odd support, while the best 
estimators of odd derivatives have even support. These estimators are "best" 
in the sense that their lowest-order error terms have a small coefficient and 
that they do not attenuate the higher frequencies as much as the alternative 
ones. A good estimator of the second derivative of z, for example, is 



{zxx}k « -jCzjfc-i - 2-2* + z*+i)> 
while a good estimator of the first derivative of z is just 

{z x }k » -(«*+! - z k ). 



(92) 



(93) 



Note that the latter, like other estimators with even support for odd deriva- 
tives, gives an estimate valid at the point midway between samples. 

This suggests that one should use staggered grids. That is, the arrays 
containing sampled values of p and q (and hence image brightness E) should 
be offset by 1/2 pixel in both * and y from those for z (see Figure 3). This 
also means that if the image is rectangular and contains n x m pixels, then 
the array of heights should be of size (n + 1) x (m + 1). Appropriate two- 
dimensional estimators for the first partial derivatives of z then are (see also 
[Horn & Schunck 81]): 

{z x }k,l « 2^(2*,/+ i ~ Zk < 1 + **+*.'+ 1 ~ **+u) 

i z y}k,l « 2^(zk+i,i ~ Zk,i + Zk+1,1+1 - z k ,i+i) 
These can be conveniently shown graphically in the form of the stencils: 



(94) 



_1_ 
2e 



-1 


+1 


-1 


+1 



and 



2e 



+ 1 


+1 


-1 


-1 



The results obtained apply to the point (k + 1/2, /+ 1/2) in the grid of discrete 
values of z; or the point (k, /) in the (offset) discrete grid of values of p and q. 
Similar schemes can be developed for the first partial derivatives of p and q 
needed in the algorithms introduced here, with the offsets now acting in the 
opposite direction. 
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Figure 3. It is convenient to have the discrete grid for p, q (and hence 
for the image E itself) offset by 1/2 pixel in x and 1/2 pixel in y from the 
grid for z. 



6.2 Discrete Estimators of the Laplacian 

We also need to obtain local averages based on discrete approximations of the 
Laplacian operators. We could simply use one of the stencils 
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or 
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The second, diagonal, form has a higher coefficient on the lowest-order error 
term than the first, edge-adjacent form, and so is usually not used by itself. 
The diagonal form is also typically not favored in iterative schemes for solving 
Poisson's equations, since it does not suppress certain high frequency compo- 
nents. We can write a stencil for a linear combination of the edge-adjacent 
and the diagonal versions in the form 



(a + l)e2 
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A judiciously chosen weighted average, namely one for which a = 1/5, is 
normally preferred, since this combination cancels the lowest-order error term. 
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If we wish to prevent the iterative scheme from "walking away" from the 
solution, however, we need to make our estimate of the Laplacian consistent 
with repeated application of our estimators for the first partial derivatives. 
That is, we want our discrete estimate of Az to be as close as possible to our 
discrete estimate of 

(*«)* + (*,),. (95) 

It is easy to see that the sum of the convolution of the discrete estimator for 
the x-derivative with itself and the convolution of the discrete estimator for 
the y-derivative with itself yields the diagonal pattern. So, while the diagonal 
pattern is usually not favored because it leads to less stable iterative schemes, 
it appears to be desirable here to avoid inconsistencies between discrete es- 
timators of the first and second partial derivatives. Experimentation with 
various linear combinations bears this out. The edge-adjacent stencil is very 
stable and permits over-relaxation (SOR) with a = 2 (see next section), but 
leads to some errors in the solution with noisefree input data. The diagonal 
form is less stable and requires a reduced value for cv, but allows the scheme to 
converge to the exact algebraic solution to problems that have exact solutions. 

The incipient instability inherent in use of the diagonal form is a reflection 
of the fact that if we think of the discrete grid as a checkerboard, then the 
"red" and the "black" squares are decoupled 30 . That is, updates of red squares 
are based only on existing values on red squares, while updates of black squares 
are based only on existing values on black squares. Equivalently, note that 
there is no change in the incremental update equations when we add a discrete 
function of the form 

Sz kl = (-1)*+' (96) 

to the current values of the height. The reason is that the estimators of the 
first derivatives and the diagonal form of the Laplacian estimator are com- 
pletely insensitive to components of this specific (high) spatial frequency 31 . 
Fortunately, the iterative update cannot inject components of this frequency 
either, so that if the average of the values of the "red" cells initially matches 
the average of the values of the "black" cells, then it will continue to do so. 
The above has not turned out to be an important issue, since the iteration 



30 The "red" and "black" squares are the cells for which the sum of the row and 
column indexes are even and odd respectively. 

It may appear that this difficulty stems from the use of staggered grids. The 
problem is even worse when aligned grids are used, however, because the discrete 
estimator of the Laplacian consistent with simple central difference estimators of 
the first partial derivatives has a support that includes only cells that are 2e away 
from the center. And this form of the Laplacian operator is known to be badly 
behaved. We find that there are four decoupled subsets of cells in this case. 
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appears to be stable with the diagonal form of the average, that is, for a = 1, 
when the natural boundary conditions are implemented with care. 

6.3 Boundary Conditions 

The boundary conditions have also to be dealt with properly to assure con- 
sistency between first- and second-order derivative estimators. In a simple 
rectangular image region, the natural boundary conditions for z could be 
implemented by simply taking the average of the two nearest values of the 
appropriate gradient component and multiplying by e to obtain an offset from 
the nearest value of z in the interior of the grid. That is, for 1 < k < n and 

1 < / < m, we could use 

Zk,0 = Zk,\ - ^(pk-i,o +Pk,o) 

Zk,m = Zk,m-1 + j:(Pk-l,m-l + Pk,m-l) 

e 2 (97) 

zo,i = zi,/ - ^(qo,l-i + qo,l) 

z n ,l = z n -i,i + -(q n -i,i-i + q n -i,l) 

on the left, right, bottom and top border of a rectangular image region (the 
corners are extrapolated diagonally from the nearest point in the interior using 
both components of the gradient). But this introduces a connection between 
the "red" and the "black" cells, and so must be in conflict with the underlying 
discrete estimators of the derivatives that are being used. 

One can do better using offsets from cells in the interior that lie in diag- 
onal direction from the ones on the boundary. That is, for 2 < k < n — 1 and 

2 < / < m — 1, we use 

Zk,o = g ( z k-i,i ~ e(Pk-i,o - qk-i,o) + Zk+i,i - e(p k ,o + qk,o)) 

Zk,m = X ( z k-l,m-l + t{pk-l,m-\ + tfit-l.m-l) + 2fc+l,m-l + e(Pk,m-l ~ tfjfc.m-l)) 

zo,i = - (21,1-1 + <Po,/-i - qo,i-i) + 21,7+1 - e(po,l + qo,l)) 

z n ,l = -(z n -i,/-i +e(p„-i,/_i +q n -i t i-i) -l-Zn-i^+i - e(p n _ ij / - q n -i,l)) 

(98) 
on the left, right, bottom and top border of a rectangular image region. The 
corners are again extrapolated diagonally from the nearest point in the interior 
using both components of the gradient. Note that, in this scheme, one point 
on each side of the corner has to be similarly interpolated, because only one 



6. Some Implementation Details 35 

of the two values needed by the above diagonal template lies in the interior 
of the region. 

If the surface gradient is not given on the image boundary, then natural 
boundary conditions must be used for p and q as well. The natural boundary 
condition is that the normal derivatives of p and q are zero. The simplest 
implementation is perhaps, for 1 < k < n — 1 and 1 < I < m — 1, 

Pk,o = Pk,l 

Pk,m-1 = Pk,m-2 

PO,l = Pi, I 

Pn-\,l =Pn-2,l 

and similarly for q (points in the corner are copied from the nearest neighbor 
diagonally in the interior of the region). It may be better to again use a 
different implementation, where the values for points on the boundary are 
computed from values at interior cells that have the same "color." That is, 
for 2 < k < n - 2 and 2 < / < m - 2, 

Pk,o = gCpt-i.i + Pk+l,l), 
1 r 

Pk,m-1 — -{Pk-l,m-2 + Pk+l,m-2), 

i (ioo) 

Po,i = ^(Pi.'-i +Pi,i+i), 

Pn-l,l = -{p n -2,l-\ +p n -2,l+l), 

and similarly for q. As before, the corner points, and one point on each side 
of the corner have to be copied diagonally, without averaging, since only one 
of the two values needed lies in the interior of the region. 

6.4 Iterative Schemes and Parallelism 

There are numerous iterative schemes for solution of large sparse sets of equa- 
tions, among them: 

Gauss-Seidel — with replacement — sequential update; 
Jacobi — without replacement — parallel update; 
Successive Over-Relaxation; 
Kazmarz relaxation; 
Line relaxation. 

Successive over-relaxation (SOR) makes an adjustment from the old value 
that is a times the correction computed from the basic equations. That is, 
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for example, 



4? +1) = 4? ) +«(4? ) -4? ) ) (ioi) 

where z y kl is the "new" value calculated by the ordinary scheme without over- 
relaxation. When a > 1, this amounts to moving further in the direction of 
the adjustment than suggested by the basic equations. This can speed up 
convergence, but also may lead to instability 32 . The Gauss-Seidel method 
typically can be sped up in this fashion by choosing a value for a close to 
two — the scheme becomes unstable for a > 2. Unfortunately the Gauss-Seidel 
method does not lend itself to parallel implementation. 

The Jacobi method is suited for parallel implementation, but successive 
over- relaxation cannot be applied directly — the scheme diverges for a > 1. 
This greatly reduces the speed of convergence. Some intuition may be gained 
into why successive over-relaxation cannot be used in this case, when it is 
noted that the neighbors of a particular cell, the ones on which the future value 
of the cell is based, are changed in the same iterative step as the cell itself. 
This does not happen if we use the Gauss-Seidel method, which accounts for 
its stability. This also suggests a modification of the Jacobi method, where 
the parallel update of the cells is divided into sequential updates of subsets 
of the cells. Imagine coloring the cells in such a way that the neighbors of a 
given cell used in computing its new value have a different color from the cell 
itself. Now it is "safe" to update all the cells of one color in parallel (for an 
analogous solution to a problem in binary image processing, see chapter 4 of 
[Horn 86]). 

Successive over-relaxation can be used with this modified Jacobi method. 
If local averages are computed using only the four edge-adjacent neighbors of 
a cell, then only two colors are needed (where the colors are assigned according 
to whether i + j is even or odd — see Figure 4). Each step of the iteration is 
carried out in two sub-steps, one for each of the cells of one color. The above 
shows that the improved convergence rates of successive over-relaxation can 
be made accessible to parallel implementations. 

When the illumination of the surface is oblique (light source away from 
the viewer), R(p, q) will tend to be locally approximately linear. This means 
that the gradient of R(p, q) will point in more or less the same direction over 
some region of the image. The effect of this is that influences on the adjust- 
ments of the estimated gradient tend to be much smaller along a direction at 
right angles to the direction "away from the light source," than they are along 
other directions. This can be seen most easily when the coordinate system 
is aligned with the direction toward a single light source in such a way that 

Conversely, if the basic method has a tendency to be unstable, then one can 
"under- relax" — that is, use a value a < 1. 
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Figure 4. The modified Jacobi method operates on subsets of cells with 
different "colors" at different times. In the simplest case, there are only 
two colors, one for the cells where the sum of the indexes is even, the other 
for the cells where the sum of the indexes is odd. 



the reflectance map has bilateral symmetry with respect to the axis q = 0. 
Then R q will be small, at least for gradients near the p-axis. In this case the 
coefficients on the diagonal of the 2x2 matrix may be very different in mag- 
nitude. This is analogous to a system of equations being much stiffer in one 
direction than another, and suggests that the convergence rate may be lower 
in this case. A possible response to this difficulty is the use of line relaxation. 

6.5 Aliasing, and How to Avoid It 



Discrete samples can represent a continuous waveform uniquely only if the con- 
tinuous waveform does not contain frequency components above the Nyquist 
rate (u>o = ir/e, where e is the spacing between samples). If a waveform is 
sampled that contains higher frequency components, these make contributions 
to the sampled result that are not distinguishable from low frequency compo- 
nents. If, for example, we have a component at frequency uq < u < 2u?o, it 
will make the same contributions as a component at frequency 2o; — u>. This 
is what is meant by aliasing. Ideally, the continuous function to be sampled 
should first be lowpass filtered. Filtering after sampling can only suppress 
desirable signal components along with aliased information. 

Numerical estimation of derivatives is weakly ill-posed. The continuous 
derivative operator multiplies the amplitude of each spatial frequency com- 
ponent by the frequency, thus suppressing low frequencies and accentuating 
higher frequencies. Any corruption of the higher frequencies is noticeable, 
particularly if most of the signal itself is concentrated at lower frequencies. 
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This means that we have to be careful how we estimate derivatives and how 
we sample the image. 

Suppose, for example, that we have an image of a certain size, but that we 
would like to run our shape-from-shading algorithm on a smaller version, per- 
haps to obtain a result in a reasonable amount of time, or to cheaply provide 
useful initial values for iteration on the finer grid. It would be quite wrong 
to simply sub-sample the original image. Simple block-averaging is better, al- 
though frequency analysis shows that the response of a block-averaging filter 
first drops to zero only at twice the Nyquist frequency. It is better to use a 
cubic spline approximation of the ideal 

sin(7rx/e) ^ 102 ^ 

(■Kx/e) 
response for filtering before sub-sampling [Rifman & McKinnon 74] [Bern- 
stein 76] [Keys 81] [Abdou k Young 82]. There is nothing specific in the above 
relating to shape-from-shading; these are considerations that apply generally 
to machine vision. 

Similar notions apply to processing of the surface itself. If we have a 
digital terrain model of a certain resolution and want to generate a lower 
resolution shaded image from it, we need to first filter and sample the digital 
terrain model. Otherwise the result will be subject to aliasing, and some 
features of the shaded image will not relate in a recognizable way to features 
of the surface. 

Finally, in creating synthetic data it is not advisable to compute the 
surface gradient on a regular discrete set of points and then use the reflectance 
map to calculate the expected brightness values. At the very least one should 
perform this computation on a grid that is much finer than the final image, 
and then compute block averages of the result to simulate the effect of finite 
sensing element areas — just as is done in computer graphics to reduce aliasing 
effects 33 . 

(This hints at an interesting problem, by the way, since the brightness 
associated with the average surface orientation of a patch is typically not quite 
equal to the average brightness of the surface, since the reflectance map is not 
linear in the gradient. This means that one has to use a reflectance map 
appropriate to the resolution one is working at — the reflectance map depends 
on the optical properties of the micro-structure of the surface, and what is 
micro-structure depends on what scale one is viewing the surface at.) 



33 0ne can obtain good synthetic data, however, with an exact algebraic solution, 
by sampling the height on a regular discrete set of points and then estimating 
the derivatives numerically, as discussed in section 5.1. This was done here to 
generate most of the examples shown in section 7. 
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6.6 Measuring the Quality of Reconstruction 

There are many ways of accessing the quality of the solution surface generated. 
Not all are useful: 

• In the case of a synthetic image obtained from a surface model, the best 
test of the output of a shape-from-shading algorithm is comparison of the 
surface orientations of the computed result with those of the underlying 
surface. One can either compute the root-mean-square deviation of the 
directions of the computed surface normals from the true surface normals, 
or just the root-mean-square difference in the gradients themselves. 

• Shading is a function of the surface gradient and thus most sensitive 
to higher spatial frequencies. Conversely, in the presence of noise and 
reconstruction errors, we expect that the lower spatial frequencies will not 
be recovered as well as the higher ones. This makes pointwise comparison 
of the heights of the computed surface with that of the original surface 
somewhat less useful, since errors in the lower spatial frequencies will 
affect this result strongly. Also, errors in height will be a function of the 
width of the region over which one has attempted to recover height from 
gradient. 

• Comparison of an "image" obtained by making brightness a function 
of height with a similar "image" obtained from the original surface is 
usually also not very useful, since such a representation is not sensitive 
to surface orientation errors, only gross errors in surface height. Also, 
people generally find such displays quite hard to interpret. 

• Oblique views of "wire-meshes" or "block-diagrams" defined by the dis- 
crete known points on the surface may be helpful to get a qualitative idea 
of surface shape, but can be misleading and are difficult to compare. If 
the shape-from-shading scheme is working anything like it is supposed 
to, the differences between the solution and the true surface are likely to 
be too small to be apparent using this mode of presentation. 

• Comparing the original image with an image obtained under the same 
lighting conditions from the solution is not useful, since the brightness 
error is reduced very quickly with most iterative schemes. Also, a "so- 
lution" can have gradient field {pkl,<lkl} that yields exactly the correct 
image when illuminated appropriately, yet it may not even be integrable. 
In fact, the "surface" may yield an arbitrary second image when illumi- 
nated from a different direction unless p and q are forced to be consistent 
(that is, unless p y = q x ) as discussed at the end of section 7.3. 

• If the underlying surface is known, shaded views of the solution and 
the original surface, produced under lighting conditions different from 
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those used to generate the input to the algorithm, are worth comparing. 
This is a useful test, that immediately shows up shortcomings of the 
solution method. It also is a graphic way of portraying the progress 
of the iteration — one that is easier to interpret than a set of numbers 
representing the state of the computation. 

• Various measures of departure from integrability may be computed. Per- 
haps most useful are comparisons of numerical estimates of (z x ,z y ) with 
(p, q). Slightly less useful is the difference (p y — q x ) of the solution, since 
the height z may still not have converged to the best fit to p and q, even 
when the gradient itself is almost integrable. 

6.7 When to Stop Iterating 

As is the case with many iterative processes, it is difficult to decide when to 
stop iterating. If we knew what the underlying surface was, we could just wait 
for the gradient of the solution to approach that of the surface. But, other 
than when we test the algorithm on synthetic images, we do not know what 
the surface is, otherwise we would probably not be using a shape-from-shading 
method in the first place! Some other tests quantities include: 

• The brightness error 

JJ(E(x,y)-R(p,q)) 2 dxdy (103) 

should be small. Unfortunately this error becomes small after just a few 
iterations, so it does not yield a useful stopping criterion. 

• The departure from smoothness 

JJ(pl+pl + ql + q 2 y )dxdy (104) 

also drops as the solution is approached, but it does not constitute a par- 
ticularly good indicator of approach to the solution. In particular, when 
one comes close to the solution, one may wish to reduce the parameters 
A, perhaps even to zero, in which case further iterations may actually 
reduce smoothness in order to better satisfy the remaining criteria. 

• One of the measures of lack of integrability 

(Py-Qx) 2 dxdy (105) 

is also not too useful, since it can at times become small,or stop changing 
significantly, even when z is still inconsistent with p and q. 

• Another measure of lack of integrability 

//((** " Pf + (*» - 9?) dx dy (106) 
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appears to be very useful, since it drops slowly and often keeps on chang- 
ing until the solution has converged. 
• One can also keep track of the rate of change of the solution with itera- 
tions 

//(i) ,+ (!)'**■ < io7) 

One should not stop until this has become small. In most cases it helps 
to continue for a while after the above measures stop changing rapidly 
since the solution often continues to adjust. 

Some of the implementation details given above may appear to be extraneous. 
However, when all of these matters are attended to, then the iterative algo- 
rithm will not "walk away" from the solution, and it will find the solution, 
to machine precision, given exact data (and assuming that boundary condi- 
tions for p and q are given, and that A' is reduced to zero as the solution is 
approached). Convergence to the exact solution will not occur when some- 
thing is amiss, such as a mismatch between the discrete estimators of the first 
derivative and the discrete estimator of the Laplacian. It is not yet clear how 
significant all of this is when one works with real image data, where there is no 
exact solution, and where the error introduced by incorrect implementation 
detail may be swamped by errors from other sources. 

7. Some Experimental Results 

The new algorithm has been applied to a number of synthetic images of simple 
shapes (such as an asymmetrical Gaussian, a sum of Gaussian blobs, and a 
sum of low frequency sinusoidal gratings) generated with a number of different 
reflectance maps (including one linear in p and q, Lambertian with oblique 
illumination, and a rotationally symmetric one). These synthetic images were 
small (usually 64 x 64 picture cells) in order to keep the computational time 
down. Typically the surface normals would be within a degree or two of the 
correct direction after a few hundred iterations. With appropriate boundary 
conditions, the computed shape would eventually be the same (to machine 
precision) as the shape used to generate the synthetic image. In each case, 
the brightness error decreased rapidly, while the integrability of the estimated 
gradient decreased much more slowly. 

7.1 Graphical Depiction of Solution Process 

For help in debugging the algorithm, and for purposes of determining a good 
schedule for adjusting the parameters \i and A', it is useful to print out the 
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(run-schedule 2) 

Lanbda: 1.0000 Mu: 8.1000 Eps: 1.0000 filpha-pq: 

Iter: Grad E: 2.09619 Bright E: 0.52570 Int-z 

Iter: 4 Grad E: 0.60296 Bright E: 0.17914 Int-z 

Iter: 8 Grad E: 0.27103 Bright E: 8.07793 Int-z 

Lanbda: 8.5000 Mu: 0.1000 Eps: 1.0000 Rlpha-pq 

Iter: 14 Grad E: 0.21348 Bright E: 0.03485 Int-z 

Iter: 20 Grad E: 8.28982 Bright E: 8.03317 Int-z 

Iter: 26 Grad E: 0.20821 Bright E: 0.03287 Int-z 

Iter: 32 Grad E: 0.20782 Bright E: 0.03268 Int-z 

Lanbda: 8.2888 Mu: 8.1800 Eps: 1.0000 Rlpha-pq: 

Iter: 40 Grad E: 0.18785 Bright E: 8.02212 Int-z 

Iter: 48 Grad E: 0.18496 Bright E: 0.02183 Int-z 

Iter: 56 Grad E: 0.18334 Bright E: 8.82174 Int-z 

Iter: 64 Grad E: 8.18226 Bright E: 0.02169 Int-z 

Lanbda: 0.1000 Mu: 0.1000 Eps: 1.0000 filpha-pq: 

Iter: 80 Grad E: 0.16684 Bright E: 0.81523 Int-z 

Iter: 96 Grad E: 8.16446 Bright E: 8.81516 Int-z 

Iter: 112 Grad E: 0.16301 Bright E: 0.01514 Int-z 



Iter: 128 Grad E: 8.16195 Bright E: 0.01513 Int-z 
Lanbda: 0.0500 Mu: 0.1008 Eps: 1.0080 Rlpha-pq: 
Iter: 160 Grad E: 8.14655 Bright E: 0.01022 Int-z 
Iter: 192 Grad E: 0.14411 Bright E: 0.01821 Int-z 
Iter: 224 Grad E: 0.14252 Bright E: 0.01021 Int-z 
Iter: 256 Grad E: 0.14139 Bright E: 0.01021 Int-z 
Lanbda: 0.0288 Mu: §.8888 Eps: 1.0000 Rlpha-pq: 
Iter: 320 Grad E: 0.12432 Bright E: 0.00573 Int-z 
Iter: 384 Grad E: 0.11074 Bright E: 0.00394 Int-z 
Iter: 448 Grad E: 0.10616 Bright E: 8.00426 Int-z 
Iter: 512 Grad E: 0.10092 Bright E: 0.88381 Int-z 
Lanbda: 8.0180 Mu: 0.0508 Eps: 1.0000 Rlpha-pq: 
Iter: 640 Grad E: 0.09924 Bright E: 0.00378 Int 
Iter: 768 Grad E: 0.09640 Bright E: 0.00376 Int-z E 
Iter: 896 Grad E: 0.09519 Bright E: 0.00374 Int-z E 
Iter:1824 Grad E: 0.09459 Bright E: 8.80373 Int-z E 
Lanbda: 0.3050 Mu: 8.8388 Eps: 1.8888 Rlpha-pq: 1 
Iter:1280 Grad E: 0.18298 Bright E: 0.80202 Int-z E: 
Iter:1536 Grad E: 0.18338 Bright E: 8.88281 Int-z E: 
It«r:1792 Grad E: 0.10336 Bright E: 0.00201 Int-z E: 
Iter:2848 Grad E: 0.18348 Bright E: 8.80201 Int-z E: 
Lanbda: 0.0028 Mu: 0.0208 Eps: 1.8880 Rlpha-pq: 1 
Iter:2384 Grad E: 0.09030 Bright E: 0.00111 Int-z E 
Iter:2560 Grad E: 8.08954 Bright E: 8.00110 Int-z E 
Iter:2816 Grad E: 8.88932 Bright E: 8.88118 Int-z E 
Iter:3872 Grad E: 0.08924 Bright E: 0.00110 Int-z E 
Lanbda: 0.0010 Mu: 8.0100 Eps: 1.0000 filpha-pq: 
Iter: 3328 Grad E: 8.88869 Bright E: 0.00086 Int-z E 
Iter: 3584 Grad E: 8.88866 Bright E: 0.00086 Int-z E 
Iter: 3848 Grad E: 8.88865 Bright E: 0.00086 Int-z E 
Iter:4896 Grad E: 8.08864 Bright E: 8.80086 Int-z E 
Lanbda: 8.0005 Mu: 0.0088 Eps: 1.0000 Rlpha-pq: 1 
Iter:4352 Grad E: 0.07698 Bright E: 0.80050 Int-z E: 
Iter:4608 Grad E: 0.07620 Bright E: 8.00049 Int-z E: 
Iter:4864 Grad E: 8.07597 Bright E: 0.00049 Int-z E: 
Iter:5120 Grad E: 0.07588 Bright E: 0.00049 Int-z E: 
Lanbda: 0.0083 Mu: 0.0078 Eps: 1.8000 Rlpha-pq: 1 
Iter:5376 Grad E: 0.06687 Bright E: 8.00030 Int-z E: 
Iter:5632 Grad E: 0.06620 Bright E: 8.00838 Int-z E: 
Iter:5888 Grad E: 8.06599 Bright E: 8.80838 Int-z E: 
Iter:6144 Grad E: 8.06590 Bright E: 0.00030 Int-z E: 
Lanbda: 0.0002 Mu: 0.0058 Eps: 1.8888 Rlpha-pq: 1 
Iter:6656 Grad E: 8.86486 Bright E: 8.00025 Int-z E 
Iter:7168 Grad E: 8.86399 Bright E: 8.88825 Int-z E 
Iter:7688 Grad E: 8.86398 Bright E: 8.88825 Int-z E 
Iter:8192 Grad E: 8.86397 Bright E: 8.88825 Int-z E 
Lanbda: 8.8801 Mu: 0.0020 Eps: 1.0000 Rlpha-pq: : 
Iter:8704 Grad E: 0.06941 Bright E: 0.00032 Int-z E 
Iter:9216 Grad E: 0.06955 Bright E: 0.88032 Int-z E 
Iter:9728 Grad E: 8.86958 Bright E: 0.88832 Int-z E 
Iter:18248 Grad E: 0.06958 Bright E: 0.00032 Int-z E 
Lanbda: 0.0000 Mu: 8.8818 Eps: 1.0000 Rlpha-pq: 1 
Iter:10752 Grad E: 0.01967 Bright E: 0.00001 Int-z E 
Iter:11264 Grad E: 0.00742 Bright E: 0.00000 Int-z E 
Iter:11776 Grad E: 0.00366 Bright E: 0.00000 Int-z E 
Iter:12288 Grad E: 0.00189 Bright E: 0.00000 Int-z E 



.7000 Rl 
0.00000 
1.38119 
0.32173 

.7000 Rl 
0.08287 
0.86189 
8.05496 
0.05077 

.7880 Rl 
0.04226 
0.03885 
8.83689 
0.03562 

.7000 Rl 
0.02934 
0.02881 
8.02736 
0.32696 

.7000 Rl 
0.02194 
0.02155 
a. 821 38 
8.82128 
7000 Rl 
0.01768 
0.01262 
0.01341 
0.01188 

,7800 Rl 
0.01167 
0.01153 
0.01147 
0.81144 

,7000 Rl 
0.01598 
8.81592 
8.01590 
0.01590 

,7088 Rl 
0.81341 
0.81337 
0.01336 
0.01336 

.7000 Rl 
0.01434 
0.01433 
0.81433 
0.01433 

,7000 Rl 
0.01111 
0.01108 
0.01107 
0.01107 

,7000 Rl 
0.00890 
0.00887 
0.00887 
8.88886 

.7888 Rl 
8 . 80861 
8.88868 
8.80860 
8.80868 
7000 Rl 
0.01012 
0.01012 
0.01012 

0.01812 
7000 Rl 
0.88186 
0.00039 
0.00017 
0.88008 



pha-z: 1.7000 Cor-pq 0.2500 Cor 



pha-z: 1.7088 Cor-pq 0.2508 Cor 



pha-z: 1.7808 Cor-pq 8.2588 Cor-: 



pha-z: 1.7000 Coi — pq 0.2500 Cor 



pha-z: 1.7000 Cor-pq 0.2508 Cor 



pha-z: 1.7888 Cor-pq 8.2500 Cor 



pha-z: 1.7888 Cor-pq 8.2500 Cor 



Int-pq E: 0.88880 
Int-pq E: 8.47459 
Int-pq E: 8.89847 



Int-pq E: 8.82512 

Int-pq E: 8.81987 

Int-pq E: 0.01918 

Int-pq E: 0.01886 



Int-pq E: 8.82146 

Int-pq E: 8.82111 

Int-pq E: 8.82891 

Int-pq E: 8.02077 



Int-pq E: 0.02160 

Int-pq E: 0.02143 

Int-pq E: 0.02136 

Int-pq E: 8.82131 



Int-pq E: 8.82062 

Int-pq E: 8.82056 

Int-pq E: 8.82052 

Int-pq E: 8.82849 



Int-pq E: 8.81893 

Int-pq E: 8.81359 

Int-pq E: 0.01483 

Int-pq E: 8.01326 



Int-pq E: 0.01318 

Int-pq E: 0.81388 

Int-pq E: 8.81382 

Int-pq E: 8.81299 
pha-z: 1 

Int-pq E: 8.81858 

Int-pq E: 8.81841 

Int-pq E: 8.81839 

Int-pq E: 8.81839 



pha-z: 1.7000 Cor-pq 0.2508 Cor 



Int-pq E: 8.81576 

Int-pq E: 8.81578 

Int-pq E: 8.81568 

Int-pq E: 8.81567 



pha-z: 1.7000 Coi — pq 0.2500 Cor- 



Int-pq E: 8.81679 

Int-pq E: 8.81679 

Int-pq E: 8.81678 

Int-pq E: 8.01678 



pha-z: 1.7000 Cor-pq 0.2500 Cor 



Int-pq E: 8.81298 

Int-pq E: 8.01285 

Int-pq E: 0.01284 

Int-pq E: 8.01283 



pha-z: 1.7000 Cor-pq 0.2500 Cor 



Int-pq E: 0. 01026 

Int-pq E: 0.01022 

Int-pq E: 0.01021 

Int-pq E: 0.01020 



pha-z: 1.7000 Cor-pq 0.2500 Cor 
Int-pq E: 0.00988' 
Int-pq E: 0.00988 
Int-pq E: 0.08988 
Int-pq E: 8.88988 



Unsnooth: 
Unsnooth: 
Unsnooth : 6 
Unsnooth: 6 
pha-z: 1.7888 Cor-pq 8.2500 Cor-: 



Unsnooth: 3 
Unsnooth: 0. 
Unsnooth: 8 



Unsnooth: 8 

Unsnooth: 8 

Unsnooth: 

Unsnooth: 



Unsnooth : 

Unsnooth: 0. 

Unsnooth: 0. 

Unsnooth: 0. 



Unsnooth: 0. 

Unsnooth: 8. 

Unsnooth: 6. 

Unsnooth: 8. 



Unsnooth: 8. 

Unsnooth: 0. 

Unsnooth: 0. 

Unsnooth: 0. 



Unsnooth: 0. 

Unsnooth : . 

Unsnooth : . 

Unsnooth: 8. 



Unsnooth: 0. 

Unsnooth: 8. 

Unsnooth: 0. 

Unsnooth: 0. 
7000 Cor-pq 0.2500 Cor- 

Unsnooth: 0. 

Unsnooth: 0. 

Unsnooth: 0. 

Unsnooth: 0. 



Unsnooth: 0. 

Unsnooth: 0. 

Unsnooth: 0. 

Unsnooth: 0. 



Unsnooth: 

Unsnooth: 

Unsnooth: 0. 

Unsnooth: 0. 



Unsnooth: 0. 

Unsnooth: 0. 

Unsnooth: 0. 

Unsnooth: 0. 



Unsnooth: 0. 

Unsnooth: 0. 

Unsnooth: 0. 

Unsnooth: 8. 



Int-pq E: 8.81162) Unsnooth 
Int-pq E: 8.81163) Unsnooth: 8. 
Int-pq E: 8.81163) Unsnooth: 8. 
(Int-pq E: 8.81163) Unsnooth: 8 
pha-z: 1.7888 Cor-pq 8.2588 Cor-: 
(Int-pq E: 8.80092) Unsnooth: 
(Int-pq E: 0.00038) Unsnooth: 8 
(Int-pq E: 8.88814) Unsnooth: 8 
(Int-pq E: 0.00007) Unsnooth: 



z 1.0000 (64 X 64) 

24685 

66932 dpq: 1.091 998 dz: 2. 173221 

19669 dpq: 0.292286 dz: 0.522765 

z 1.0888 (64 X 64) 

86248 dpq: 0.044399 dz: 0.070357 

85529 dpq: 8.886952 dz: 8. 02401 9 

05508 dpq:0. 002019 dz:0. 016154 

85518 dpq: 8. 881 283 dz: 8. 81 2275 

z 1.8888 (64 X 64) 

86969 dpq: 8. 88251 3 dz: 8. 81 1044 

07009 dpq: 0.000822 dz: 8. 007575 

07019 dpq: 8.888681 dz: 8. 006001 

07022 dpq: 8.888485 dz: 8. 885886 

z 1.0000 (64 X 64) 

08152 dpq:8. 888535 dz:8. 834017 

08167 dpq:B. 888324 dz:8. 082881 

08170 dpq:0. 000252 dz:B. 082305 

08170 dpq:0. 0O02O4 dz:0. 001920 

z 1.0000 (64 X 64) 

09233 dpq:0. 000219 dz:8. 881581 

89246 dpq:B. BBB142 dz:B. 881074 

09253 dpq:0. 000107 dz:0. 880839 

09257 dpq:0. 000884 dz:0. 000673 

z 1.0000 (64 X 64) 

10434 dpq:0. 000101 dz:0. 000698 

11383 dpq:0. 000586 dz:0. 001413 

11413 dpq:0. 010930 dz:0. 004026 

11442 dpq: 0. 000078 dz: 0.088526 

z 1.8888 (64 X 64) 

11467 dpq:B. 800068 dz:0. 000432 

11514 dpq: 8. 880033 dz: 0.000251 

11537 dpq:0. 800021 dz:0. 000164 

11551 dpq:0. 800014 dz:0. 000111 

z 1.8888 (64 X 64) 

11612 dpq: 8. 883722 dz: 0.000633 

11608 dpq: 0. 0031 55 dz: 0.008624 

11608 dpq: 0.003287 dz: 0.000635 

11687 dpq: 8. 8031 41 dz: 0.080629 

z i.eeee (64 x 64) 

12242 dpq: 8. 887526 dz: 8. 881322 
12259 dpq:8. 887359 dz:8. 801289 
12264 dpq:8. 887298 dz:8. 881277 
12267 dpq: 8. 807276 dz: 0.001273 
z 1.008B (64 X 64) 

12387 dpq:0. 011513 dz:8. 802398 

12388 dpq:8. 811497 dz:8. 882395 
12388 dpq: 8. 81 1493 dz: 8. 882394 
12388 dpq:B. 811492 dz:8. 882394 
z 1.00BB (64 X 64) 

12855 dpq:8. 887626 dz:8. 881396 
12875 dpq:8. 887546 dz:8. 881378 
12882 dpq:B. 887589 dz:8. 881371 
12885 dpq:B. 887494 dz:8. 881368 
z 1.0000 (64 X 64) 
13258 dpq:0. 005249 dz:B. 888933 
13278 dpq:B. 885163 dz:8. 888917 
13285 dpq: 8. 8851 19 dz: 8. 888989 
13288 dpq: 0.005899 dz: 8.888986 
z 1.0000 (64 X 64) 
13378 dpq:0. 0O5282 dz:0. 000968 
13381 dpq: 8. 885278 dz: 8. 888966 

13381 dpq: 0.O05268 dz: 0.000965 

13382 dpq: 0.085268 dz: 8. 888965 
z 1.8880 (64 X 64) 
13186 dpq: 8. 887721 dz: 8. 881 538 
13182 dpq:0. 007740 dz:0. 001535 
13181 dpq:0. 007744 dz:8. 881536 
.13188 dpq: 8. 887744 dz: 8. 881536 
z 1.8888 (64 X 64) 
.16893 dpq:0. 000062 dz:B. B8B185 
.16351 dpq:8. 000813 dz:B. 888062 
.16446 dpq:0. 000005 dz:0. 000025 

16488 dpq:0. 000003 dz:0. 000011 



Figure 5. Diagnostic trace of various error measures. This sequence of 
results corresponds to the reconstruction of the sharp-edged crater shape 
shown in Figure 7. This kind of presentation is important, but must be 
supplemented by some graphic depiction of the evolving solution surface. 
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diagnostic measurements discussed in sections 6.6 and 6.7. But it is hard 
to tell exactly what is going on just by looking at a large table of numbers 
such as that shown in Figure 5. It is important to also provide some graphic 
depiction of the evolving shape as it is computed. To make shaded images 
of the reconstructed surface useful, however, they must be illuminated from 
a direction different from the direction of illumination used for the original 
input image 34 . Shown in Figure 6, is such a sequence of shaded images gener- 
ated during the reconstruction of the surface of a polyhedral object, starting 
from a random field of surface orientations. Here the image provided to the 
algorithm 35 corresponded to illumination form the Northwest, while illumi- 
nation from the Northeast was used to display the reconstruction. Note how 
the edges become sharper as A', controlling the contribution of the penalty 
term for departure from smoothness, is made smaller and smaller. This ex- 
ample also illustrates the algorithm's ability to deal with surfaces that have 
discontinuities in surface orientation. 

Because of the interest in application to astrogeology, a crater-like shape 
was also reconstructed, as shown in Figure 7. In this case, the algorithm 
rapidly found a shape that was generally correct, except for flaws in places on 
the rim of the crater in the Northwest and Northeast. These are areas where 
there is little contrast between the inside and the outside of the crater in the 
input image 36 . It took the algorithm a considerable number of additional 
iterations to determine the correct continuation of the shape computed in 
other image areas. 

7.2 Emergent Global Organization 

Often progress towards the correct solution is not as uneventful. Frequently 
small internally consistent solution patches will establish themselves, with 
discontinuities in surface orientation where these patches adjoin. Also, coni- 
cal singularities form that tend to move along the boundaries between such 
regions as the iterative solution progresses. Conversely, boundaries between 
solution patches often form along curves connecting conical singularities that 
form earlier. After a large enough number of iterations, patches of local orga- 
nization tend to coalesce and lead to emergent global organization. This can 

34 The test illumination should be quite different from the illumination used to 

generate the original image — preferrably lying in a direction that differs from the 

original source direction by as much as n/2 
35 The input image is not shown, but is just like the last image in the sequence 

shown, except that left and right are reversed. 
36 Again, the input image is not shown, but is like the last image in the sequence 

shown, except that left and right are reversed. 
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Figure 6. Reconstruction of a portion of a truncated hexahedron from 
a shaded image. The dihedral angle between each pair of the three outer 
surfaces is w/2. This example illustrates the algorithm's ability to deal 
with surfaces that have large discontinuities in surface orientation. 
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Figure 7. Reconstruction of a crater-like shape. Points on the 
the Northeast and the Southwest correspond to places in the input image 
where there is least contrast between the inside and the outside, since the 
direction of the incident illumination is parallel to the rim there. 
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be observed best when A' is smaller than what it would normally be set to for 
rapid convergence. In Figures 8 and 9, for example, are shown a sequence of 
shapes leading finally to a spherical cap on a planar surface. Within some re- 
gions, solution surface patches quickly establish themselves that individually 
provide good matches to corresponding parts of the input image. The bor- 
ders between these internally consistent regions provide error contributions 
that the algorithm slowly reduces by moving the boundaries and incremen- 
tally changing the shapes within each of the regions. Too rapid a reduction 
of A' can remove the incentive to reduce the creases and kinks and to freeze 
the solution in a state where some unnecessary discontinuities remain. If, for 
example, A' were to be set to zero with a "solution" consisting of a spherical 
cap with an inner disk inverted, as in the right hand image of the middle 
row of Figure 9, there would be no incentive to further reduce the length of 
the circular discontinuity, and the smooth solution for this part of the image 
would not be found. 

The algorithm was also applied to impossible shaded images. Suppose, 
for example, that we are dealing with a Lambertian surface illuminated by 
a source near the viewer and that there is a dark smudge in the middle of a 
large planar region facing us (which appears brightly lit). It turns out that 
there is no surface with continuous first derivatives that could give rise to a 
shaded image with a simply connected, bounded dark region in the middle 
of a bright region [Szeliski & Horn 89]. In Figure 10 we see what happens 
when the algorithm attempts to find a solution. Patches grow within which 
the solution is consistent with the image, but which yield discontinuities at 
boundaries between patches. Conical singularities sit astride these boundaries. 
For all random initial conditions tried, the algorithm eventually eliminates all 
but one of these conical singularities. The computed surface is in fact a 
"solution," if one is willing to allow such singularities. 

The graphical method of presenting the progress of the iterative solutions 
illustrated above was very helpful in debugging the program and in determin- 
ing reasonable schedules for reduction of the parameters A' and fi. Shown in 
Figure 1 1 are some examples of what happens when things go wrong. In the 
top row are shown instabilities arising in the solution for the crater-like shape, 
near the points where there is low contrast between the inside and the outside 
of the crater — that is, where there is no local evidence for curvature. These 
instabilities can be suppressed by reducing A' more slowly. In the middle row 
are shown patterns resulting from various programming errors. Finally, in the 
bottom row is shown the propagation of an instability from a free boundary 
when A' is set to zero. It appears that the process is not stable when there 
are neither boundary conditions for the height nor for the gradient. 
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Figure 8. Emergent global organization of local nonlinear iterative pro- 
cess. Internally consistent solutions arise in certain image patches with 
discontinuities at the borders between regions. The boundaries between 
these patches move, and the solutions within the patches adjust in order 
to reduce the sum of the error terms. 
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Figure 9. Emergent global organization of local nonlinear iterative pro- 
cess. Neighboring patches coalesce, as conical singularities are absorbed 
by coalescing with other singularities, or by being pushed towards a con- 
tour where the surface orientation is discontinuous. The final solution is 
a spherical cap resting on a plane. 
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Figure 10. What happens when the algorithm is confronted with an 
impossible shaded image? The input image here (not shown) is a circularly 
symmetric dark smudge in a uniformly bright surround. The light source 
is assumed to be near the viewer. The algorithm finds a "solution" with 
a single conical singularity. 
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Figure 11. Graphical depiction of instabilities and the effects of pro- 
gramming errors. In the top row are shown instabilities resulting from too 
rapid reduction of the penalty term for departure from smoothness. The 
middle row shows the results of various programming errors. The bottom 
row shows waves of instability propagating inwards from a free boundary. 
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Figure 12. New one-step shape-from-shading algorithm. See text! 
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In the past, shape-from- shading algorithms have often been "tested" by 
verifying that the computed gradient field actually generates something close 
to the given input image. To show just how dangerous this is, consider Fig- 
ure 12, which demonstrates a new non-iterative method for recovering a "sur- 
face" given a shaded image. In Figure 12(a), is the input to the algorithm, 
while Figure 12(c) is what the gradient field that is constructed by this al- 
gorithm looks like when illuminated in the same way as the original surface. 
Now Figure 12(b) shows what the original surface looks like when illuminated 
from another direction. As a test, we should check that the computed gradi- 
ent field looks the same under these illuminating conditions. But behold, it 
does not! In Figure 12(d) we see what we get when we use the other illumi- 
nating condition. The "trick" here is that the problem of shape from shading 
is heavily underconstrained if we are only computing a gradient field and not 
enforcing integrability. There are many solutions and we can, in fact, impose 
additional constraints. The underlying gradient field here was computed by 
solving the photometric stereo equations [Woodham 78, 79, 80, 89] for the two 
images in Figures 12(c) and (d) under the two assumed lighting conditions. 

The new algorithm has also been applied to synthetic images generated 
from more complicated surfaces such as digital terrain models (DTMs) used 
earlier in research on interpretation of satellite images of hilly terrain [Horn Sz 
Bachman 78] [Sjoberg &: Horn 83] and in automatic generation of shaded over- 
lays for topographic maps [Horn 81]. These synthetic images were somewhat 
larger (the one used for Figure 1 is of size 231 x 178, for example). In this 
case, the simple algorithm, presented in section 5.3, using a regularizing term 
would often get trapped in a local minimum of the error function after a small 
number of iterations, while the modified algorithm presented in section 5.6, 
exploiting the linearization of the reflectance map, was able to proceed to a 
solution to machine precision after a few thousand iterations. Most of the 
surface normals typically were already within a degree or so of the correct 
direction after a few hundred iterations. 

The closeness of approach to the true solution depends on several of the 
implementation details discussed earlier. In particular, it was helpful to use 
the old values of p and q for the reference point in the linearization of R(p, q), 
rather than any of the other choices suggested earlier. Also, it helps to use 
the diagonal averaging scheme in the iteration for height rather than the one 
based on edge-adjacent neighbors. 

7.3 Rating the Difficulty of Shape-from-Shading Problems 

Experiments with synthetic shaded images suggests that certain shape-from- 
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shading problems are relatively easy, while others are quite difficult. First 
of all, basso-relievo 37 surfaces (those with only low slopes) are easy to deal 
with (see also section 2.6) in comparison with alto-relievo surfaces (those with 
steep slopes). The digital terrain model used for the experiment illustrated 
in Figure 1 falls in the latter category, since the sides of the glacial cirque are 
steep and the individual gullies steeper still. 

Typically the brightness of a surface patch increases the more it is turned 
towards the light source. If it is turned too far, however, it becomes so steep 
that its brightness once again decreases. There is a qualitative difference 
between shape-from-shading problems where none of the surface patches are 
turned that far, and those where some surface patches are so steep as to have 
reduced brightness. In the latter case, there appears to be a sort of two- 
way ambiguity locally about whether a patch is dark because it has not been 
turned enough to face the light source or whether it has been turned too far. 
This ensures that simplistic schemes will get trapped in local minima where 
patches of the solution have quite the wrong orientation. Similarly, the more 
sophisticated scheme described here takes many more iterations to unkink the 
resulting creases. 

The transition between the two situations depends on where the light 
source is. The difficulty is reduced when the illumination is oblique (see also 
section 2.6). Conversely, the problem is more severe when the light source is 
at the viewer, in which case brightness decreases with slope independent of 
the direction of the surface gradient. This explains why the algorithm took 
longer to find the solution in the case of the spherical cap (Figure 8) since 
it was illuminated by a source near the viewer. It was more straightforward 
to find the solutions for the truncated hexahedron and the crater-like surface 
(Figures 6 and 7), both of which where illuminated obliquely. The above 
dichotomy is related to another factor: problems where the relevant part of 
the reflectance map is nearly linear in gradient are considerably easier to deal 
with than those in which the reflectance map displays strong curvatures of 
iso-brightness contours. 

Smooth surfaces, particularly when convex, can be recovered easily. Sur- 
faces with rapid undulations and wrinkles, such as the digital terrain model 
surface (Figure 1) are harder. Discontinuities in surface orientation are even 
more difficult to deal with. Note that, with the exception of the digital terrain 
model, all of the examples given here involve surfaces that have some curves 
along which the surface orientation is not continuous. The spherical cap, for 
example, lies on a planar surface, with a discontinuity in surface orientation 



For more regarding the terms basso-relievo, mezzo-relievo and alto-relievo, see 
[Koenderink & van Doom 80]. 
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where it touches the plane. 

Problems where boundary conditions are not available, and where there 
are no occluding boundaries or singular points, are ill-posed. Not too surpris- 
ingly these tend to lead to instabilities in the algorithm, particularly when 
one attempts to reduce the penalty term for departure from smoothness. In 
these cases instabilities can be damped out to some extent by enforcing the 
image irradiance equation on the boundary by iterative adjustment of the 
gradient computed from the discrete approximation of the natural boundary 
conditions for p and q. But results have not been promising enough to discuss 
here in more detail. 

The number of iterations to converge to a good solution appears to grow 
almost quadratically with image size (number of rows or columns). This is 
because some effects have to "diffuse" across the image. This means that the 
total amount of computation grows almost with the fourth power of the (lin- 
ear) image size. It is well known that ordinary iterative schemes for solving 
elliptic partial differential equations quickly damp out higher spatial frequency 
errors, while low frequency components are removed very slowly. The way to 
deal with this problem is to use computation on coarser grids to reduce the 
low spatial frequency components of the error. This is the classic multigrid 
approach [Brandt 77, 80, 84]. It is clear that a true multigrid implementation 
(as opposed to a simple pyramid scheme) 38 would be required to pursue this 
approach further on larger images. This is mostly to cut down on the com- 
putational effort, but can also be expected to reduce even further the chance 
of getting caught in a local minimum of the error function. Implementation, 
however, is not trivial, since the equations are nonlinear, and because there 
are boundary conditions. Both of these factors complicate matters, and it is 
known that poor implementation can greatly reduce the favorable convergence 
rate of the basic multigrid scheme [Brandt 77, 80, 84]. 

Alternatively, one can apply so-called direct methods for solving Poisson's 
equations [Simchony, Chellappa & Shao 89]. 

8. Conclusion 

The original approach to the general shape-from-shading problem requires 
numerical solution of the characteristic strip equations that arise from the 
first-order nonlinear partial differential equation that relates image irradiance 
to scene radiance [Horn 70, 75]. Variational approaches to the problem instead 



38 A naive approach has one solve the equations on a coarse grid first, with the 
results used as initial conditions for a finer grid solution after interpolation. True 
multigrid methods are more complex, but also have much better properties. 
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minimize the sum of the brightness error and a penalty term such as a measure 
of departure from smoothness. These yield second-order partial differential 
equations whose discrete approximation on a regular grid can be conveniently 
solved by classic iterative techniques from numerical analysis. Several of these 
methods, however, compute surface orientation, not height, and do not ensure 
that the resulting gradient field is integrable [Ikeuchi & Horn 81] [Brooks & 
Horn 85]. One thus has, as a second step, to find a surface whose gradient 
comes closest to the estimated gradient field in a least-squares sense (see 
[Ikeuchi 84], chapter 11 in [Horn 86], and [Horn & Brooks 86]). 

The two steps can be combined, and the accuracy of the estimated surface 
shape improved considerably, by alternately taking one step of the iteration for 
recovering surface orientation from brightness, and one step of the iteration 
that recovers the surface that best fits the current estimate of the surface 
gradient. This idea can be formalized by setting up a variational problem 
involving both the surface height above a reference plane and the first partial 
derivatives thereof. The resulting set of three coupled Euler equations can be 
discretized and solved much as the two coupled equations are in the simpler 
methods that only recover surface orientation. 

Such an iterative scheme for recovering shape from shading has been 
developed and implemented. The new scheme recovers height and gradient at 
the same time. Linearization of the reflectance map about the local average 
surface orientation greatly improves the performance of the new algorithm 
and could be used to improve the performance of existing iterative shape- 
from-shading algorithms. The new algorithm has been successfully applied to 
complex wrinkled surfaces; even surfaces with discontinuities in the gradient. 
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